© E. Vitral and G. A. Mamon 2021

1. Introduction

When the Laser Interferometer Gravitational-Wave Observatory (LIGO) first detected gravitational waves coming from a stellar-mass black hole merger (Abbott et al. 2016) and then the Event Horizon Telescope (EHT) released the first image of the supermassive black hole (SMBH) in M 87 (Event Horizon Telescope Collaboration 2019), astronomers obtained the most compelling evidence about the existence of those intriguing and particular objects. Yet, black holes (BHs) have been treated as more than a theoretical object for a considerable amount of time, starting in 1939, when Oppenheimer & Snyder (1939) proposed them to be the final step of the life of massive stars (≳10 M_⊙) after their final gravitational collapse into stellar-mass black holes, with masses in between ∼3 M_⊙ (Thompson et al. 2020) and ≈52 M_⊙ (Woosley 2017), but also when Hoyle & Fowler (1963) identified the then recently discovered quasars (Schmidt 1963) as SMBHs. This latter class of BHs, which reside in the centers of massive galaxies, are responsible for extremely luminous sources in the Universe, such as quasars and active galactic nuclei (AGN), sometimes unleashing powerful jets of relativistic matter, along with outflows that have a profound impact on star formation and galaxy evolution (e.g., Croton et al. 2006; Hopkins et al. 2006). In addition, some black holes may have formed during the early moments of the Universe (e.g., Zel’dovich & Novikov 1966; Hawking 1971), and these primordial BHs may constitute an import mass fraction of black holes below the SMBH mass.

Since there is no theoretical constraint to the mass of a black hole, it would be reasonable to believe that intermediate-mass black holes (IMBHs) could exist, filling the considerable gap between stellar-mass black holes and SMBHs (i.e., with masses between 100 and 10⁵ M_⊙). Furthermore, SMBHs are understood to grow by mergers, where the first seeds are stellar-mass BHs (Madau & Rees 2001) or metal-free primordial gas clouds (Loeb & Rasio 1994), so IMBHs may be a transitory stage in the growth of BHs. However, there is currently little evidence for IMBHs (see reviews by Volonteri 2010; Greene et al. 2020), with some important candidates highlighted (e.g., Kaaret et al. 2001; Chilingarian et al. 2018 in dwarf galaxies, and recently Lin et al. 2020 in a globular cluster) and one gravitational wave confirmation (The LIGO Scientific Collaboration & the Virgo Collaboration 2020). Furthermore, IMBHs could help to explain many enigmas in astrophysics, such as filling up part of the dark matter mass budget (e.g., Haehnelt & Rees 1993; Loeb & Rasio 1994) or providing massive seeds for high redshift quasars, whose high masses at such early times represent a challenge to current theories (Haiman 2013). Therefore, great efforts have been undertaken to detect IMBHs to better understand their origin and evolution.

Globular clusters (GCs) appear to be a unique laboratory to test the existence of IMBHs. These quasi-spherical star clusters are known to have old stellar populations, indicating that they formed at early epochs. Their high stellar number densities provide an excellent environment to increase stellar interactions that could give birth to compact objects. More precisely, contrary to galaxies, the rates of stellar encounters in the inner parts of GCs containing half their stellar mass are sufficiently high to statistically affect the orbits of their stars by two-body relaxation (Chandrasekhar 1942). Moreover, after several relaxation times, the interplay between the negative and positive heat capacities of the inner core and outer envelope causes transfers of energy. This in turn leads to the gravothermal catastrophe, where the core collapses resulting in a steep inner density profile, while the envelope expands. Roughly one-fifth of GCs are believed to have suffered such core-collapse (Djorgovski & King 1986).

Several scenarios have been proposed for the existence of IMBHs in GCs. One is the direct collapse of population III stars (Madau & Rees 2001), but the link of population III stars with GCs is not clear. Another is the accretion of residual gas on stellar-mass BHs formed in the first generation of stars (Leigh et al. 2013), but the availability of the gas is unclear as the first massive stars will blow it out of the GC by supernova explosions. Stellar mergers are a popular mechanism for IMBH formation. Portegies Zwart & McMillan (2002) proposed a runaway path to IMBH formation in GCs, where an initially massive star suffers multiple physical collisions with other stars during the first few Myr of the GC, before they have time to explode as supernovae or simply lose mass. During these collisions, the most massive star will lose linear momentum, ending up at the bottom of the gravitational potential well. At the same time, its mass will grow during the successive stellar mergers, to the point that it will end up as a BH, possibly reaching 0.1% of the GC stellar mass. Miller & Hamilton (2002) proposed a slower process for IMBH formation, where dynamical friction (Chandrasekhar 1943) causes the most massive stellar remnant BHs to sink to the center of the gravitational well over Gyr. Thus a ≳50 M_⊙ stellar remnant BH, sufficiently massive to avoid being ejected from the GC by dynamical interactions, would grow in mass through mergers with these other massive BHs as well as other typically massive stars, reaching a mass of 1000 M_⊙ over the Hubble time, which they argued generates IMBHs in some ten percent of GCs. Finally, Giersz et al. (2015) proposed that hard binaries containing stellar-mass BHs merge with other stars and binaries, which can be a fast or slow process.

These models present, however, drawbacks: The short relaxation time needed in the Portegies Zwart & McMillan (2002) scenario usually requires primordial mass segregation in order not to eliminate too many GCs candidates, while the assumption by Miller & Hamilton (2002) of BH seeds above ≈50 M_⊙ is not expected as the massive progenitors are fully exploded in pair-instability supernovae (e.g., Woosley 2017).

Unfortunately, attempts to detect IMBHs have been somewhat inconclusive: Dynamical modeling is still dependent on the assumptions concerning the confusion between the IMBH and a central subcluster of stellar remnants (e.g., den Brok et al. 2014; Mann et al. 2019; Zocchi et al. 2019). Furthermore, these analyses usually rely on too few stars inside the sphere of influence of the IMBH, which can lead to false detections (Aros et al. 2020). Besides, searches for signs of accretion indicate no strong evidence for > 1000 M_⊙ black holes in galactic GCs (e.g., Tremou et al. 2018).

In this paper, we analyze the core-collapsed Milky Way globular cluster NGC 6397, and search for kinematic imprints of a central IMBH or subcluster of unresolved objects (hereafter, CUO). We use very precise and deep Hubble Space Telescope (HST) proper motions (PMs) from Bellini et al. (2014) and line-of-sight (LOS) velocities from the Multi Unit Spectroscopic Explorer (MUSE) instrument on the Very Large Telescope, to which we added PMs from the Gaia Data Release 2 (hereafter, Gaia or Gaia DR2).

NGC 6397 has been broadly investigated in previous studies. The PMs lead to plane-of-sky (POS) velocities with equal radial and tangential dispersions, leading to the appearance of isotropic orbits (Heyl et al. 2012; Watkins et al. 2015a with HST for the inner orbits, and Jindal et al. 2019 with Gaia for the outer orbits). But this projected velocity isotropy can hide three-dimensional (3D) variations of the anisotropy of the 3D velocity ellipsoid. Kamann et al. (2016) used the JAM mass-orbit modeling code (Cappellari 2008) to fit the LOS velocities obtained with the MUSE spectrograph, and found that a 600 ± 200 M_⊙ black hole best fits their data, agreeing with the radio emission upper limit of 610 M_⊙ for the same cluster provided in Tremou et al. (2018). On the other hand, numerical N-body modeling by Baumgardt (2017) strongly excluded the possibility of an IMBH in NGC 6397, not to say that much of the non-luminous mass measured in this cluster could actually be in the form of unresolved white dwarfs and low-mass stars (Heggie & Hut 1996).

Many of the mass modeling studies performed with this cluster assumed or constrained its surface density parameters to values estimated long ago (e.g., Trager et al. 1995), that may suffer from problems such as radial incompleteness, which is better addressed by missions such as HST and Gaia. This lack of accuracy on the surface density can strongly impact the dynamical analysis of NGC 6397.

In the present study, we performed state-of-the-art mass-orbit modeling to analyze the presence of an IMBH or a CUO in the center of NGC 6397 and better measure the radial variation of its velocity anisotropy. Among the different popular mass/orbit modeling (e.g., Mamon et al. 2013; Watkins et al. 2013; Read & Steger 2017; Vasiliev 2019a; see Chap. 5 of Courteau et al. 2014 for a review), we used an extension of the MAMPOSST: Modeling Anisotropy and Mass Profiles of Observed Spherical Systems (Mamon et al. 2013) Bayesian code, which now takes into account both LOS velocities and PMs (MAMPOSST-PM, Mamon & Vitral, in prep.). MAMPOSST-PM was found to perform extremely well in a code challenge using mock data for dwarf spheroidal galaxies, with many similarities to GCs (Read et al. 2020).

We profit from the versatility of MAMPOSST-PM to fit not only cases with a single population of stars, but also with two mass populations, thus allowing for mass segregation. We also provide new fits to the surface density parameters of this cluster by jointly modeling HST and Gaia data.

2. Method

2.1. MAMPOSSt-PM

We summarize here the main aspects of MAMPOSST-PM (details in Mamon & Vitral, in prep.). MAMPOSST-PM is a Bayesian code that fits parametrized forms of the radial profiles of mass and velocity anisotropy to the distribution of observed tracers (here, GC stars) in two-, three- or four-dimensional projected phase space (PPS): Projected distances to the GC center (hereafter projected radii), combined with LOS velocities and/or PMs. The mass profile is the sum of observed tracers, as well as a possible central BH, and additional dark components. By working on discrete stars, MAMPOSST-PM can probe the inner regions better than methods where velocity or PM data are binned in radial intervals (e.g., van der Marel & Anderson 2010), and an analysis of dark matter on mock dwarf spheroidals has shown that the conclusions can depend on how the data are binned (Richardson & Fairbairn 2014).

The probability density that a tracer at projected radius R has a velocity vector v is

$$\begin{array}{c}\hfill p(\mathit{v}|R)=\frac{g(R,\mathit{v})}{\int g(R,\mathit{v})\mathrm{d}\mathit{v}}=\frac{g(R,\mathit{v})}{\mathrm{\Sigma }(R)}\end{array}$$(1)

where Σ(R) is the surface density (SD), while g(R, v) is the density in PPS, such that

$$\begin{array}{cc}\hfill {\displaystyle \int \int }2\pi Rg(R,\mathit{v})\mathrm{d}R\mathrm{d}\mathit{v}& =\mathrm{\Delta }{N}_{\mathrm{p}}\hfill \\ \hfill & =N_{\mathrm{p}}(R_{\max })-N_{\mathrm{p}}(R_{\min }),\hfill \end{array}$$(2)

with $N_{\text{p}}(R)={\int }_0^{R}2\pi R^{\prime }\mathrm{\Sigma }(R^{\prime })\text{d}{R}^{\prime }$ being the projected number of points (e.g., stars in a star cluster) in a cylinder of radius R. The likelihood is then calculated as

$$\begin{array}{c}\hfill \mathcal{L}={\displaystyle \prod _{i}p({\mathit{v}}_i|R_i),}\end{array}$$(3)

where the i indices represent individual points (stars).

The observed stars tend to be located in in the GC. But, inevitably, some observed stars will be interlopers (here field stars, FS) mostly foreground in our case, but possibly some background too. MAMPOSST-PM splits the PPS distribution into separate GC and FS components. Equation (1) becomes

$$\begin{array}{cc}\hfill p(\mathit{v}|R)& =\frac{g_{\mathrm{GC}}(R,\mathit{v})+g_{\mathrm{FS}}(R,\mathit{v})}{{\mathrm{\Sigma }}_{\mathrm{GC}}(R)+{\mathrm{\Sigma }}_{\mathrm{FS}}(R)}\hfill \\ \hfill & =\frac{g_{\mathrm{GC}}(R,\mathit{v})+{\mathrm{\Sigma }}_{\mathrm{FS}}{f}_{\mathrm{LOS}}(v_{\mathrm{LOS}})f_{\mathrm{PM}}(\mathit{PM})/{\eta }^2}{{\mathrm{\Sigma }}_{\mathrm{GC}}(R)+{\mathrm{\Sigma }}_{\mathrm{FS}}},\hfill \end{array}$$(4)

where PM is the vector of PMs (μ_α, *, μ_δ), f_LOS(v_LOS) and f_PM(PM) are the respective distribution functions of field star LOS velocities and PM vectors¹ while η = 4.7405 D is the conversion of PM in mas yr⁻¹ to POS velocity in km s⁻¹ given the distance D to the system, in kpc. The second equality of Eq. (4) assumes that the FS surface density is independent of position.

The GC contribution to the PPS density is the mean local GC velocity distribution function h, averaged along the LOS (Mamon et al. 2013; Read et al. 2020):

$$\begin{array}{c}\hfill g_{\mathrm{GC}}(R,\mathit{v})=2{\displaystyle {\int }_{\mathrm{R}}^{\infty }h(\mathit{v}|R,r)\nu (r)\frac{r}{\sqrt{r^2-R^2}}\mathrm{d}r,}\end{array}$$(5)

where r is the 3D distance to the system center, while ν(r) is the 3D tracer density (here stellar number density) profile.

MAMPOSST-PM assumes that the local velocity ellipsoid, in other words the local 3D velocity distribution function (VDF), is separable along the three spherical coordinates (r, θ, ϕ):

$$\begin{array}{c}\hfill h(\mathit{v}|R,r)=h(v_r|R,r)h(v_{\theta }|R,r)h(v_{\varphi }|R,r).\end{array}$$(6)

MAMPOSST-PM further assumes that the three one-dimensional local VDFs are Gaussian:

$$\begin{array}{c}\hfill h(v_i|R,r)=\frac{1}{\sqrt{2\pi {\sigma }_i^2}}exp(-\frac{v_i^2}{2{\sigma }_i^2})·\end{array}$$(7)

This Gaussian assumption for the local VDFs is much better than the very popular Gaussian assumption for the LOS-integrated VDF, p(v|R) because velocity anisotropy affects the shape of p(v|R), in particular that of p(v_LOS|R) (Merritt 1987). By measuring the shape of the PPS, hence of p(v|R), MAMPOSST-PM gets a good handle on velocity anisotropy (see below), and therefore on the mass profile because of the mass-anisotropy degeneracy (Binney & Mamon 1982). MAMPOSST-PM considers the measurement errors by adding them in quadrature to the velocity dispersion for the GC component.

The velocity anisotropy profile β at radius r is defined as

$$\begin{array}{c}\hfill {\displaystyle \beta (r)=1-\frac{{\sigma }_{\theta }^2(r)+{\sigma }_{\varphi }^2(r)}{2{\sigma }_r^2(r)},}\end{array}$$(8)

where θ and ϕ are the tangential components of the coordinate system and r is the radial component, while ${\sigma }_i^2$ stands for the velocity dispersion of the component i of the coordinate system.

When considering spherical symmetry in velocity space, the variances of the composite Gaussian VDF h(v|(R, r) can be written as

$$\begin{array}{cc}& {\sigma }_{\mathrm{LOS}}^2(R,r)=[1-\beta (r){\left(\frac{R}{r}\right)}^2]{\sigma }_{\mathrm{r}}^2(r),\hfill \end{array}$$(9a)

$$\begin{array}{cc}& {\sigma }_{{\mathrm{POS}}_{\mathrm{R}}}^2(R,r)=[1-\beta (r)+\beta (r){\left(\frac{R}{r}\right)}^2]{\sigma }_{\mathrm{r}}^2(r),\hfill \end{array}$$(9b)

$$\begin{array}{cc}& {\sigma }_{\mathrm{POSt}}^2(R,r)=[1-\beta (r)]{\sigma }_{\mathrm{r}}^2(r),\hfill \end{array}$$(9c)

where Eq. (9a) is from Binney & Mamon (1982), while Eqs. (9b) and (9c) are from Strigari et al. (2007). In Eqs. (9a)–(9c), the radial velocity variance ${\sigma }_{\text{r}}^2(r)$ is obtained by solving the spherical stationary Jeans equation with no streaming motions (in particular, rotation):

$$\begin{array}{c}\hfill \frac{\mathrm{d}\left(\nu {\sigma }_{\mathrm{r}}^2\right)}{\mathrm{d}r}+2\frac{\beta (r)}{r}\nu (r){\sigma }_{\mathrm{r}}^2(r)=-\nu (r)\frac{GM(r)}{r^2},\end{array}$$(10)

where M(r) is the total mass profile.

The distribution of interlopers in PPS is straightforward: The spatial density is assumed to be uniform (hence Σ_FS does not depend on R in Eq. (4)). The distribution of LOS velocities, f_LOS(v_LOS) in Eq. (4), is assumed Gaussian. The two-dimensional distribution of PMs, f_PM(PM) in Eq. (4), is found to have wider wings than a Gaussian (see Sect. 5.2.3 and Appendix B.1). The POS velocity errors are not added in quadrature to the velocity dispersions because the PM distribution of interlopers has wider tails (Appendix B.2). Instead, MAMPOSST-PM uses an approximation (Mamon & Vitral, in prep.) for the convolution of PM errors with the PM dispersions.

2.2. Dark component

MAMPOSST-PM allows for a diffuse dark component with a large choice of analytical density profiles. This dark component could represent dark matter or alternatively unseen stars.

Dark matter (DM) may dominate the outskirts of GCs, and perhaps also deeper inside. Contrary to distant GCs (Ibata et al. 2013), the presence of substantial DM in the outskirts of NGC 6397 is unlikely because this GC appears to have an orbit around the Milky Way that takes it from an apocenter of roughly 6.6 kpc to a pericenter of 2.9 kpc in 88 Myr (Gaia Collaboration 2018a), so it should have suffered from dozens of previous episodes of tidal stripping, leaving little DM (Mashchenko & Sills 2005), if there was any in the first place. In particular, the orbital parameters of Gaia Collaboration (2018a), Z = −0.48 ± 0.01 kpc and W = −127.9 ± 3 km s⁻¹, indicate that NGC 6397 just passed through the Milky Way disk only ≈4 Myr ago².

One may ask whether there is substantial DM in the inner regions of GCs like NGC 6397. Using N-body simulations of isolated ultra-compact dwarf galaxies, which appear to be larger analogs of GCs, but with relaxation times longer than the age of the Universe, Baumgardt & Mieske (2008) found that the inner DM component is heated by dynamical friction from the stars, causing a shallower inner DM distribution than that of the stars. Shin et al. (2013) traced the dynamical evolution of NGC 6397 with Fokker-Planck methods, discarding those particles that reach the GC tidal radius. They found that the GC could have contained up to one-quarter of its initial mass in the form of DM and match the present-day surface density and LOS velocity dispersion profiles.

On the other hand, there may be unseen matter in the core of the GC that is not dark matter per se, but composed of unseen stars, for example an inner nuclear cluster of faint stars, possibly white dwarfs, neutron stars, or stellar-mass black holes (Zocchi et al. 2019; Mann et al. 2019) and binary stars (Mann et al. 2019). We thus performed MAMPOSST-PM runs including such a CUO instead of (or in addition to) an IMBH.

2.3. Velocity anisotropy profile

With LOS data, MAMPOSST is able (at least partially) to lift the degeneracy first pointed out by Binney & Mamon (1982) between the radial profiles of mass and velocity anisotropy (Mamon et al. 2013). With the three components of the velocity vector, MAMPOSST-PM is even much more able to disentangle mass and anisotropy profiles (Read et al. 2020; Mamon & Vitral in prep.). With the data used in this paper, we have thus a unique chance of constraining the mass profile of our sources, which is essential to restrain the estimates on the IMBH.

While GCs are often modeled with isotropic velocities, we performed many runs of MAMPOSST-PM with freedom in the anisotropy profile. Indeed, a central IMBH may modify the orbital shapes. Furthermore, the outer orbits of isolated GCs are often thought to be quasi radial (Takahashi 1995). Moreover, the frequent passages of NGC 6397 across the Galactic disk (Sect. 1) could alter the orbital shapes of GC stars.

The anisotropic runs of MAMPOSST-PM used the generalization (hereafter gOM) of the Osipkov-Merritt model (Osipkov 1979; Merritt 1985) for the velocity anisotropy profile:

$$\begin{array}{c}\hfill {\displaystyle {\beta }_{\mathrm{gOM}}(r)={\beta }_0+({\beta }_{\infty }-{\beta }_0)\frac{r^2}{r^2+r_{\beta }^2},}\end{array}$$(11)

where r_β is the anisotropy radius, which can be fixed as the scale radius of the luminous tracer by MAMPOSST-PM³.

3. Global structure and distance of NGC 6397

Our mass-orbit modeling of NGC 6397 (Sect. 2 below) assumes a distance, the knowledge of the GC center, spherical symmetry, and the absence of rotation. We investigate these assumptions for NGC 6397.

3.1. Spherical symmetry

NGC 6397 appears to be close to spherical symmetry. Its elongation on the sky is 0.07 (Harris 2010).

3.2. Center

Our analysis assumes that the IMBH is located at the GC center, so our choice of center is critical. We have considered two centers: (RA, Dec) = ($265\stackrel{°}{.}175375$, $-53\stackrel{°}{.}674333$) (Goldsbury et al. 2010, using HST data) and (RA, Dec) = ($265\stackrel{°}{.}1697$, $-53\stackrel{°}{.}6773$) (Gaia Collaboration 2018a, using Gaia data), both in epoch J2000.

We selected the center with the highest local stellar counts. It is dangerous to use our HST star counts for this because the positions are relative to the center, hence need to be rescaled assuming a given center (see Sect. 4.1). We used instead the Gaia star counts, which benefit from absolute positional calibration. Figure 1 shows a Gaussian-smoothed map of Gaia star counts in the inner 100″ × 100″ (half the HST field of view) of NGC 6397, where we overplotted the centers obtained by Goldsbury et al. (2010) and Gaia Collaboration (2018a).

Clearly, the center of Gaia Collaboration (2018a) is less well aligned with the density map than the center of Goldsbury et al. (2010). We therefore select the Goldsbury et al. (2010) center at ($265\stackrel{°}{.}1754$, $-53\stackrel{°}{.}6743$) in J2000 equatorial coordinates.

3.3. Distance

The adopted distance is important because the mass profile at a given angular radius (e.g., r in arcmin) deduced from the Jeans equation of local dynamical equilibrium (Eq. (10)) varies roughly as $M\propto r{\sigma }_v^2$, thus as distance D for given LOS velocities and as D³ for given PMs. Table 1 shows a list of distance estimates for NGC 6397. Recent estimates are bimodal around 2.39 kpc and 2.64 kpc. The CMD-based distance estimates favor the larger distance, the kinematics favor the smaller distance, while the parallax distances point to small (HST) or large (Gaia) distances.

We estimated a kinematical distance by equating the LOS velocity dispersion of stars measured with VLT/MUSE (see Sect. 4.3, below) with the HST PM dispersions of the same stars (see Sect. 4.1). To avoid Milky Way field stars, we used the cleaned MUSE and HST samples, for which we had 692 matches, among which 445 with separations smaller than $0\stackrel{″}{.}1$⁴. For these 445 stars, we measured a LOS velocity dispersion of 4.91 ± 0.16 km s⁻¹, and an HST PM dispersion of 0.439 ± 0.015 mas yr⁻¹ in the RA direction and 0.442 ± 0.015 mas yr⁻¹ along the Dec direction (where the uncertainties are taken as the values divided by $\sqrt{2(n-1)}$). This yields a kinematic distance of 4.91/[c (0.439 + 0.442)/2] = 2.35 ± 0.10 kpc, where c = 4.7405 is the POS velocity of a star of PM = 1 mas yr⁻¹ located at D = 1 kpc.

We adopted the lower distance of 2.39 kpc (i.e., a distance modulus of 11.89), for three reasons: (1) We trust more the HST parallax than the Gaia-DR2 parallax, given that the former is based on a much longer baseline; (2) Our study of NGC 6397 is based on kinematics, and thus the larger distance would lead to abnormally high POS velocity dispersions compared to the LOS velocity dispersions. (3) It is consistent with our kinematic estimate of the distance. We do not adopt our estimated kinematic distance of 2.35 ± 0.10 kpc because it is based on the perfect equality of LOS and POS velocity dispersions, which supposes velocity isotropy, which in turn is not certain. For our adopted distance of 2.39 kpc, 1 arcmin subtends 0.7 pc.

4. Data

4.1. HST data

The HST data were kindly provided by A. Bellini, who measured PMs for over 1.3 million stars in 22 GCs, including NGC 6397 (Bellini et al. 2014). The data for NGC 6397 has a 202 arcsec square field of view (Wide Field Camera) and reached down to less than 1 arcsec from the Goldsbury et al. (2010) center (see the right panel of Fig. 8). This minimum projected radius to the center is smaller than the BH radius of influence $r_{\text{BH}}~G{M}_{\text{BH}}/{\sigma }_0^2\simeq 0.11(M_{\text{BH}}/600M_{\odot })$ pc (using the LOS velocity dispersion of Sect. 3.3), which corresponds to 9 (M_BH/600 M_⊙) arcsec, given our adopted distance of 2.39 kpc.

This dataset was provided in a particular master frame shape (for details, see Table 29 from Bellini et al. 2014 as well as Anderson et al. 2008), with both GC center and PM mean shifted to zero. The first step in our analysis was to convert the positions and PMs to the absolute frame.

4.1.1. HST absolute positions

We applied the Rodrigues (1840) rotation formula to shift the relative positions back to their original center to translate the GC stars to their true positions on the sky. We used the center of Goldsbury et al. (2010), which was the one considered by Bellini et al. (2014). We then rotated the subset with respect to its true center, so that the stars originally parallel to the dataset’s increasing x axis remained parallel to the right ascension increasing direction. We then verified our method by matching the stars in sky position with Gaia (see Sect. 4.4 below).

4.1.2. HST absolute proper motions

The HST PMs were measured relative to the bulk PM of the GC. We corrected the relative PMs of Bellini et al. (2014) with their provided PM corrections. We just added Cols. 4 and 5 to Cols. 31 and 32 from Table 29 of Bellini et al. (2014), respectively. We then converted the relative PMs to absolute PMs by computing the bulk PM of NGC 6397 using the stellar PMs provided by Gaia DR2 as explained in Sect. 4.2 below.

The small field of view of HST and the few pointed observations do not allow the observation of sufficiently numerous background quasars to obtain an absolute calibration of HST PMs. On the other hand, the Gaia reference frame obtained with more than half a million quasars provides a median positional uncertainty of 0.12 mas for G <  18 stars (Gaia Collaboration 2018b) and therefore allows us, by combining its accuracy with HST’s precision, to know NGC 6397 PMs with unprecedented accuracy. We will compare the PMs of stars measured both by HST and Gaia in Sect. 4.4.

Our HST data had 13 593 stars. The PM precision varies across the field of view because of the different time baselines (principally 1.9 and 5.6 years) in part due to the smaller size of the older HST cameras. Figure 2 shows that the PMs are much more accurate in a large square portion of the field of view, including the western part and extending to the central region, where the baseline is 5.6 years. At magnitude F606W = 17, the one-dimensional PM precision is 0.02 mas yr⁻¹ in the more precise portion of field of view and 0.08 mas yr⁻¹ outside.

4.2. Gaia data

Gaia DR2 presented an overall astrometric coverage of stellar velocities, positions and magnitudes of more than 10⁹ stars in the Milky Way and beyond. Gaia measured PMs in over 40 000 stars of magnitude G <  17 within a 1° cone around NGC 6397. These PMs have a typical precision of 0.17 mas yr⁻¹ at G = 17. Since the GaiaG and HST F606W magnitudes differ by less than 0.1, one sees that Gaia PMs are less precise than those of HST by factors of two to eight depending on the region of the HST field of view. The Gaia data was also used to infer the number density profile out to large distances from the GC center.

4.3. VLT/MUSE data

We complemented the PMs using LOS velocities that Husser et al. (2016) acquired with the MUSE spectrograph on the VLT. A mosaic of 5 × 5 MUSE pointings led to an effective square field of view of 5 arcmin on the side. The bulk LOS velocity is ⟨v_LOS⟩ = 17.84 ± 0.07 km s⁻¹ (Husser et al. 2016). In a companion article, Kamann et al. (2016) assigned membership probabilities to the stars according to their positions in the space of LOS velocity and metallicity, compared to the predictions for the field stars from the Besançon model of the Milky Way (Robin et al. 2003). The data, kindly provided by S. Kamann, contained 7130 LOS velocities, as well as the membership probabilities.

4.4. Consistency between datasets

4.4.1. Positional accuracy

We used Tool for OPerations on Catalogues And Tables (TOPCAT, Taylor 2005) to check the match positions between catalogs, which we did in our own Python routines that filtered the original datasets (Sect. 5 below). Performing symmetric matches among each pair of datasets, with a maximum allowed separation of 1 arcsec, we obtained median separations of $0\stackrel{″}{.}38$ for the 4455 stars in common between Gaia and HST, as well as $0\stackrel{″}{.}01$ for the 4440 stars in common between MUSE and HST.

4.4.2. Proper motions

Many of the 4455 stars in common between HST and Gaia have very high Gaia PM uncertainties (we saw in Sect. 4.2 that HST PM uncertainties were much smaller, but a few have large values). Limiting to 1179 stars in common with both Gaia and HST PM uncertainties below 0.4 mas yr⁻¹, we find that the HST PMs in the (RA, Dec) frame are (0.006 ± 0.812, 0.029 ± 1.239) mas yr⁻¹ above those from Gaia, where the “errors” represent the standard deviation. The uncertainties on the means are $\sqrt{1179}=34$ times lower: (0.024, 0.036) mas yr⁻¹. The mean shifts in PMs are (0.3, 0.8) times the uncertainties on the mean, and thus not significant. This possible offset in PMs is of no concern as long as we only use HST data, since we analyze them relative to the bulk PM of the GC. Even in MAMPOSST-PM runs with the combined HST+Gaia dataset, this shift in PMs is not statistically significant and appears too small to affect the results.

5. Data cleaning

We now describe how we selected stars from each dataset for the mass modeling runs with MAMPOSST-PM. We also used more liberal criteria to select stars in our estimates of the surface density profile (see Sect. 7.1).

For each dataset, we discarded all stars with PM errors above half of 0.394 mas yr⁻¹, corresponding (for our adopted distance) to the one-dimensional PM dispersion of NGC 6397 measured by Baumgardt et al. (2019) for the innermost 2000 Gaia stars, whose mean projected radius was $99\stackrel{″}{.}25$. The PM error is computed as the semi-major axis of the error ellipse (Eq. (B.2) of Lindegren et al. 2018):

$$\begin{array}{cc}& {ϵ}_{\mu }=\sqrt{\frac{1}{2}(C_{33}+C_{44})+\frac{1}{2}\sqrt{{(C_{44}-C_{33})}^2+4C_{34}^2}}\hfill \end{array}$$(12)

$$\begin{array}{cc}& C_{33}={ϵ}_{{\mu }_{\alpha \ast }}^2,\hfill \end{array}$$(13)

$$\begin{array}{cc}& C_{34}={ϵ}_{{\mu }_{\alpha \ast }}{ϵ}_{{\mu }_{\delta }}\rho ,\hfill \end{array}$$(14)

$$\begin{array}{cc}& C_{44}={ϵ}_{{\mu }_{\delta }}^2,\hfill \end{array}$$(15)

where ϵ denotes the error or uncertainty and where ρ is the correlation coefficient between μ_α* and μ_δ⁵. Our HST data does not provide ρ. We thus selected stars with ϵ_μ <  0.197 mas yr⁻¹, or equivalently for which the LOS velocity error satisfies ϵ_vLOS <  0.197 c D = 2.23 km s⁻¹, for D = 2.39 kpc and where c = 4.7405 (see Sect. 3.3).

Figure 3 shows the PM errors versus magnitude for HST (blue) and Gaia (red), as well as the equivalent PM error for the LOS velocity errors of Gaia (purple) and MUSE (green). The figure indicates that our maximum allowed PM errors are reached at magnitudes m_F606W = 19.7 for the worst HST data (those with the shortest baseline, see Fig. 2), and for virtually all magnitudes for the other HST stars. Gaia stars reach the maximum allowed PM error at typically G = 17.5. The equivalent LOS velocity error limit is reached at magnitude m_F606W = 16.8 for MUSE, but only at G = 13.4 for Gaia. The Gaia radial velocities are thus of little use for our modeling, as too few of them have sufficiently precise values.

5.1. HST data cleaning

After selecting the stars with low PM errors, we cleaned our HST data in three ways: we discarded stars with (1) PMs far from the bulk PM of the GC; (2) lying off the color-magnitude diagram; (3) associated with X-ray binaries.

5.1.1. HST proper motion filtering

The higher surface number density of stars in the inner regions of NGC 6397 allows us to distinguish GC stars with field stars in PM space, as shown in Fig. 4. Moreover, some high velocity stars could in principle be caused by very tight (separations smaller than ∼0.1 AU) GC binaries.

We made a very liberal cut to select GC stars in PM space, using a circle of radius 6 mas yr⁻¹ (green circle in Fig. 4), corresponding to 15 times the PM dispersion measured by Baumgardt et al. (2019). The PM center of the GC was set at the mean of the values in Table 2. Our cut in PM space also corresponds to a velocity dispersion of 75 km s⁻¹, which is over 13 times the highest LOS velocity dispersion measured by Kamann et al. (2016). We chose such a liberal cut to ensure that we would not miss any high velocity GC members because otherwise our modeling would underestimate the mass. This left us with 9149 stars among the 9624 with low PM errors.

The filtering of field stars in PM space is not fully reliable because the cloud of field stars in PM space (upper part of Fig. 4) may extend into the (smaller) GC cloud (and past it). We thus proceed to another filtering in the color-magnitude diagram (CMD).

5.1.2. HST color-magnitude filtering

The left panel of Fig. 5 shows the locations in the CMD of the stars that survived the PM error cut and the filtering in PM space. While most stars follow a tight relation, a non-negligible fraction are outliers. We took a conservative cut of the stars on the CMD using kernel density estimation as displayed in the right panel of Fig. 5⁶. This graph displays the 1, 2 and 3σ contours of the kernel density estimation, drawn as black lines. We selected stars inside the 2σ region because the 3σ region appears too wide to rule out binaries, while the 1σ region was too conservative, extracting too few stars for our analysis.

The CMD filtering not only removes field stars whose PMs coincide by chance with those of GC stars, but also removes GC members that are unresolved binary stars and lie in the edges of the main-sequence, as well as particular Blue Stragglers in the 15 <  F606W <  16.5 range, which are believed to be the result of past mergers of GC members (Leonard 1989). Removing binaries and stars who have gone through mergers is important because their kinematics are dominated by two-body interactions, while our modeling assumes that stellar motions are dominated by the global gravitational potential of the GC. As discussed by Bianchini et al. (2016), three types of binaries need to be considered.

1. Resolved (i.e., wide) binaries will produce their own PMs that can be confused with the parallax. But, following Bianchini et al. (2016), the PMs of such resolved binaries, of order of a/T where a is the semi-major axis of the binary and T is the time baseline, are negligible in comparison to the GC velocity dispersion. Indeed, to be resolved at the distance of NGC 6397, the binary star should be separated by at least $0\stackrel{″}{.}1$, which corresponds to 24 AU at the distance of NGC 6397; the most massive binaries, with mass below 2 M_⊙, will have a period of 90 years, thus a/T = 1.26 km s⁻¹, while less massive binaries will have longer periods, hence lower a/T.

2. Nearly resolved binaries will not be deblended, and their distorted image (in comparison to the PSF) will lead to less precise astrometry, which will be flagged.

3. The orbits of unresolved binaries, considered as single entities, will be those of test particles in the gravitational potential. But their higher mass should lead them to have lower velocity dispersions than ordinary stars, in particular in the dense inner regions of GCs where the two-body relaxation time is sufficiently short for mass segregation.

Our CMD filtering left us with 7259 stars among the 9149 surviving the previous filters.

5.1.3. Removal of X-ray binaries

Bahramian et al. (2020) detected 194 X-ray sources within 200 arcsec from the center of NGC 6397, some of which could potentially be background AGN. The remaining sources are thought to be X-ray binaries, and as such will have motions perturbed by their invisible compact companion. We therefore deleted the 50 HST stars whose positions coincided within 1 arcsec with the “centroid” position of an X-ray source.

Bianchini et al. (2016) found that the effect of unresolved binaries on GC dispersion profiles is only important for high initial binary fractions (e.g., ∼50%), which can induct a difference of 0.1 − 0.3 km s⁻¹ in the velocity dispersion profile, mainly in the GC inner regions, where the binary fraction should be highest. For lower initial binary fractions, Bianchini et al. (2016) found that unresolved binaries do not significantly affect the PM dispersion profile, but only the kinematics error budget. Therefore, the binary fractions calculated by Davis et al. (2008) and Milone et al. (2012a) for NGC 6397 are clearly insufficient (i.e., ≲5% and ≲7%, respectively) to require a special treatment.

5.1.4. HST final numbers

After these cuts, we are left with 7209 stars from the HST observations (among the original 13 593). The combination of the CMD and PM filtering along with the removal of X-ray binaries effectively cleans the data of most of the binaries that would affect our modeling, only leaving binaries of different luminosities that are at intermediate separations (≈0.15 to 5 AU, causing peculiar motions between 1 and 15 times the GC velocity dispersion). We checked that using a less liberal cut of 7.5σ instead of 15σ affects very little our results.

5.2. Gaia data cleaning

We followed similar steps in cleaning the Gaia data as we did for the HST data.

5.2.1. Quality flags

We filtered the Gaia stars with ϵ_μ <  0.197 mas yr⁻¹ using two data quality flags proposed by Lindegren et al. (2018). First, we only kept stars whose astrometric solution presented a sufficiently low goodness of fit:

$$\begin{array}{c}\hfill \sqrt{\frac{{\chi }^2}{N-5}}<1.2\mathrm{Max}\{1,exp[-0.2(G-19.5)]\},\end{array}$$(16)

where N is the number of points (epochs) in the astrometric fit of a given star and 5 is the number of free parameters of the astrometric fit (2 for the position, 1 for the parallax and two for the PM). Equation (16) gives a sharper HR diagram, removing artifacts such as double stars, calibration problems, and astrometric effects from binaries. It is more optimized than the (astrometric_excess_noise < 1) criterion, used in Baumgardt et al. (2019) and Vasiliev (2019b), especially for brighter stars (G ≲ 15), according to Lindegren et al. (2018).

Second, we only kept stars with good photometry.

$$\begin{array}{c}\hfill 1.0+0.015{(G_{\mathrm{BP}}-G_{\mathrm{RP}})}^2<E<1.3+0.06{(G_{\mathrm{BP}}-G_{\mathrm{RP}})}^2,\end{array}$$(17)

where E = (I_BP + I_RP)/I_G is the flux excess factor. Equation (17) performs an additional filter in the HR diagram, removing stars with considerable photometric errors in the BP and RP photometry, affecting mainly faint sources in crowded areas. This poor photometry broadens the CMD, leading to more confusion with field stars.

Those variables correspond to the following quantities in the Gaia DR2 archive:

• χ²: astrometric_chi2_al

• ν′: astrometric_n_good_obs_al

• E: phot_bp_rp_excess_factor

• G_BP − G_RP: bp_rp

5.2.2. Maximum projected radius

NGC 6397 likely suffers from tidal heating every time it passes through the Milky Way’s disk, in particular during its last passage less than 4 Myr ago (Sect. 2.2). The tidal forces felt by GC stars during passages through the disk will produce velocity impulses that are effective in perturbing the least bound orbits, which typically are those of the stars in the outer envelope of the GC. The Gaia data can trace the effects of such tidal disturbances on the GC kinematics.

We searched for anomalous kinematics in NGC 6397 using Gaia DR2 out to a maximum projected radius of 1° with respect to NGC 6397’s center. The top panel of Fig. 6 indicates that the mean POS velocities in the radial and tangential directions differ beyond 8′⁷. We also checked the concordance of the radial and tangential components of the radial profiles of POS velocity dispersion (or equivalently PM dispersion). The bottom panel of Fig. 6 shows excellent agreement between the two components of the velocity dispersion from 2′ to 20′, suggesting that the velocity ellipsoid is nearly isotropic in this range of projected radii. We conservatively adopted a maximum projected radius of 8′ as set by the divergence of the mean velocity profiles beyond that radius.

5.2.3. Gaia proper motion filtering and the bulk proper motion of NGC 6397

As for HST, we filtered Gaia stars in PM space to later filter them in CMD space. Since Gaia data extends to much greater projected radii from the GC center, thus to lower GC surface densities, the GC stands out less prominently from the field stars (FS) in PM space. We therefore first estimated the bulk PM of the GC and we assigned a first-order probability of membership using a GC+FS mixture model.

We were tempted to assign two-dimensional (2D) Gaussian distributions for both GC stars and interlopers. However, the FS PM distribution has wider tails than a Gaussian. This means that stars on the other side of the GC, relative to the center of the field star component (i.e., its bulk motion) in PM space are more likely to be field stars than assumed by the Gaussian model. We found that the PM-modulus “surface density” profile (the velocity analog of the surface density profile) is well fit by a Pearson type VII distribution (Pearson 1916), as explained in detail in Appendix B.1. This distribution relies on two free parameters, a scale radius a and an outer slope γ, and can be written as:

$$\begin{array}{c}\hfill f_{\mu }(\mathit{\mu })=-\frac{\gamma +2}{2\pi a^2}{[1+{\left(\frac{\mu }{a}\right)}^2]}^{\gamma /2},\end{array}$$(18)

where μ = (μ_α, *, μ_δ) and

$$\begin{array}{c}\hfill {\mu }_i=\sqrt{{({\mu }_{\alpha ,\ast i}-\overline{{\mu }_{\alpha ,\ast i}})}^2+{({\mu }_{\delta ,i}-\overline{{\mu }_{\delta ,i}})}^2},\end{array}$$(19)

where the suffix i stands for the component analyzed, which in the case of Eq. (18) is the interlopers (i.e., field stars, hereafter FS). The reader can verify that, indeed, ∫f_μ(μ) dμ = 1, with dμ = 2 π μ dμ.

With the respective Gaussian and Pearson VII distributions of PMs of GC stars and FS, we performed a joint fit to the two-dimensional distribution of PMs, which provided us with a precise bulk PM of the GC. For this, we considered the convolved expressions of the PM distributions of both GC and interlopers with the errors provided by the Gaia archive (ϵ_μα, *, ϵ_μδ and ρ_{μα, *μδ}). When passing onto polar coordinates, the uncertainty propagation of Eq. (19) produces

$$\begin{array}{cc}\hfill {ϵ}_{\mu ,i}^2=& {\left(\frac{{\mu }_{\alpha ,\ast i}-\overline{{\mu }_{\alpha ,\ast i}}}{{\mu }_i}\right)}^2{ϵ}_{{\mu }_{\alpha ,\ast i}}^2+{\left(\frac{{\mu }_{\delta ,i}-\overline{{\mu }_{\delta ,i}}}{{\mu }_i}\right)}^2{ϵ}_{{\mu }_{\delta ,i}}^2\hfill \\ \hfill & +2\frac{({\mu }_{\alpha ,\ast i}-\overline{{\mu }_{\alpha ,\ast i}})({\mu }_{\delta ,i}-\overline{{\mu }_{\delta ,i}})}{{\mu }_i^2}{ϵ}_{{\mu }_{\alpha ,\ast }{\mu }_{\delta ,i}},\hfill \end{array}$$(20)

where ϵ_{μα, *μδ} = ϵ_μα, * ϵ_μδ ρ_{μα, *μδ}. The convolution with Gaussian errors was straightforward in the case of GC stars since their PM distribution was also modeled as a Gaussian, and thus we just added the errors to the dispersions in quadrature:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{GC},\mathrm{new}}^2={\sigma }_{\mathrm{GC}}^2+{ϵ}_{\mu ,\mathrm{G}C}^2.\end{array}$$(21)

However, the convolution of the field star distribution with Gaussian errors cannot be reduced to an analytic function; numerical evaluation of the convolution integrals for each star would dramatically increase the calculation time. We therefore used the analytical approximation for the ratio of convolved to raw probability distribution functions of PM moduli, (which is also incorporated in MAMPOSST-PM), as briefly described in Appendix B.2 (details are given in Mamon & Vitral, in prep.). This allowed us to perform our mixture model fit to the PM data, using Markov chain Monte Carlo (MCMC)⁸ to estimate bulk motions of both the GC and the field stars and assign probabilities of GC membership for each star.

Table 2 displays our estimates of μ_α, * and μ_δ, which show a good agreement with the literature values for the GCs mean PMs. In the same table, we display the overall average of the literature plus our estimates, with its respective uncertainty ϵ, calculated as ϵ² = ⟨ϵ_i⟩² + σ², where ϵ_i stands for the uncertainties on the estimated values, and σ stands for their standard deviation. This overall average and uncertainties were also later used as priors for the MAMPOSST-PM mass modeling analysis.

We only kept stars with GC membership probability higher then 0.9 according to our mixture model. This corresponds to a 1.43 mas yr⁻¹ cut in the distribution of PMs of the GC, thus a 3.6σ cut.

5.2.4. Gaia color-magnitude filtering

We CMD-filtered the Gaia data roughly following the same KDE method as we used for the HST data. The only differences were that (1) we used the equivalent filters G, B, and R, and (2) we selected the 3σ region of the KDE, instead of 2σ, since the latter appeared to be too conservative a cut for Gaia. Figure 7 displays the final stars after this filter, along with the previously filtered subsets.

5.2.5. Removal of X-ray binaries

X-ray binaries are also in the Gaia sample. We therefore removed the five X-ray binaries that were matched to Gaia stars.

5.2.6. Gaia final numbers and comparison to HST

The cleaned Gaia sample contains 1905 stars. Figure 8 shows that the spatial coverage of Gaia DR2 is indeed poorer than HST when taking into account only stars from the clean samples. While the cleaned HST sample contains 221 stars within the 9″ IMBH radius of influence (for an IMBH with a mass of 600 M_⊙ as measured by Kamann et al. 2016, see Sect. 4.1), the cleaned Gaia presents none, which is clearly insufficient to provide reliable results on the presence of an IMBH.

5.3. MUSE data cleaning

We filtered the MUSE sample by removing stars whose probability of being a member of the GC, according to their position in velocity – metallicity space (Kamann et al. 2016, kindly provided by S. Kamann), was less than 0.9. This step left us with 6595 stars among the original 7130. We then removed stars with LOS velocity errors greater than 2.232 km s⁻¹ (half of the GC velocity dispersion), as well as stars that did not match the HST stars in a symmetric ≤1 arcsec match. We were left with 532 stars, of which 4 were previously identified as X-ray binaries (see Sect. 5.1.3), yielding a final MUSE sample containing 528 stars, thus adding LOS velocities to ∼7% of the HST sample.

Our Gaussian prior on v_LOS used the mean provided by Husser et al. (2016): ⟨v_LOS⟩ = 17.84 km s⁻¹. For the uncertainty on ⟨v_LOS⟩, we allowed a much wider dispersion of 2.5 km s⁻¹ instead of the 0.07 km s⁻¹ of Husser et al. (2016), given the relatively wide range of bulk LOS velocities reported in other studies (Milone et al. 2006; Lind et al. 2008; Carretta et al. 2009a; Lovisi et al. 2012).

5.4. Merging of the different datasets

Although we ran MAMPOSST-PM on HST and Gaia individually for data homogeneity, we preferred to combine the data from HST, MUSE and Gaia to probe the mass and velocity anisotropy parameters with more accuracy. We started with the HST filtered subset, which accounted for most of the data we would use, and then we merged it to the other datasets following the steps below:

1. We restricted the Gaia stars to those with G magnitudes within the limits of F606W magnitudes from the HST subset, (i.e., 16.11 and 22.14), given the quasi equivalence between these two filters. This step ensured that mass segregation effects would be the same for all data subsets.

2. We removed Gaia stars that were symmetrically matched to HST stars to better than 1 arcsec (see Sect. 4.4.1) since HST PMs presented smaller errors (Fig. 3).

3. We incorporated LOS velocities from MUSE, according to the approach described in Sect. 5.3.

After these steps, we were left with 8255 stars: 7209 of which had PMs from HST, with 583 of those presenting LOS velocities from MUSE, as well as 1046 additional stars with PMs from Gaia DR2.

6. Streaming motions in NGC 6397

6.1. Rotation

The presence of rotation in quasi-spherical systems makes the kinematical modeling more difficult. Any rotation will be interpreted by codes neglecting it, such as MAMPOSST-PM, as disordered motions, and should lead to different mass profiles and deduce velocity anisotropies that are more tangential than in reality.

Fortunately, NGC 6397 does not appear to have significant amount of rotation. Indeed, while Gebhardt et al. (1995) found that the integrated light of the core of NGC 6397 showed rotation with a projected amplitude of 2 km s⁻¹, Kamann et al. (2016) concluded with individual stars from VLT/MUSE observations (see Sect. 4.3 below) that rotation contributes negligibly to the second velocity moment in this GC. Vasiliev (2019c) reported a typical systematic uncertainty in Gaia DR2 PMs of up to ∼0.02 mas yr⁻¹ at any radius and considered rotation to be confirmed only when its peak amplitude exceeded ∼0.05 mas yr⁻¹, which was not the case for his NGC 6397 measurements. On the other hand, Bianchini et al. (2018) found POS rotation in NGC 6397 with Gaia DR2 at a 2σ level for this GC, but with v/σ = 0.03 only, the smallest value of the 51 GCs they analyzed. This is consistent with the quasi-null mean POSt velocity profile seen in Fig. 6. Finally, Sollima et al. (2019) detected a rotation of 0.48 km s⁻¹ from a 3D analysis, which is less than 10% of its velocity dispersion. We therefore neglect rotation in NGC 6397.

6.2. Radial motions

Similarly to rotation, any radial streaming motions will be interpreted by codes such as MAMPOSST-PM as radial dispersion. However, as seen in the top panel of Fig. 6, the Gaia DR2 PMs do not show strong signs of radial motions.

7. Practical considerations

We now discuss the practical implementation of MAMPOSST-PM to our cleaned sample of the stellar motions inside NGC 6397.

7.1. Surface density

7.1.1. Basic approach

Equation (4) requires the knowledge of the GC surface density (SD) profile, Σ_GC(R). More importantly, Eq. (5) requires the knowledge of the 3D number density ν(r) when integrating the local VDF along the LOS. While MAMPOSST-PM has a mode where it jointly fits the parameters of ν(r), M(r) and β(r) to the distribution of stars in projected phase space (see Eq. (11) of Mamon et al. 2013, which is different from our Eq. (4)), this requires the data to have constant completeness as a function of projected radius. If the astrometric measures of PMs are not independent of projected radius, one first needs to estimate (and deproject) the SD profile based on a wider dataset than that with PM values to obtain priors on the parameters of ν(r). Such an analysis had been performed for mass-orbit modeling of galaxy clusters with spectroscopic information whose completeness depended on projected radius (e.g., Mamon et al. 2019).

The PM data are likely to be incomplete at small projected radii because of the increased crowding of stars as one moves inward, especially for Gaia, whose mirrors are smaller than that of HST. We noticed that Gaia data with PMs traced less well the inner cusp of the SD than the full Gaia data (requesting only positions and G magnitude). Figure 9 shows that the Gaia PM subsample with G <  17 is indeed incomplete in the inner regions relative to the full G <  17 Gaia sample since the distribution of projected radii of Gaia stars from the GC center is shallower for stars with precise PMs than for all stars (including those without PM measurements). A Kolmogorov-Smirnov test indicates that there is only 10⁻⁵ probability that the distributions of projected radii of the subsample with precise PMs and that of the full sample (including stars without PM measurements and those with imprecise ones) arise from the same parent population. Another difficulty is to measure the outer envelope given confusion with field stars. This is important because MAMPOSST-PM computes outward integrals.

We therefore needed to provide MAMPOSST-PM with a precomputed SD profile. In MAMPOSST-PM, this profile had to be a simple analytical one- or two-parameter model (where free parameters are scale and possibly shape). More precisely, MAMPOSST-PM assumes Gaussian priors on the precomputed SD profile parameters.

7.1.2. Choice of model and main parameters

The choice of a good model for SD is crucial because the mass of a possible IMBH is linked to the inner slope of the SD profile (van der Marel & Anderson 2010). NGC 6397 has long been known to have a cuspy (steep) inner SD profile (Auriere 1982; Lauzeral et al. 1992; Lugger et al. 1995; Noyola & Gebhardt 2006; Kamann et al. 2016), which prompted Djorgovski & King (1986) to classify it as core-collapsed. Unfortunately, no simple analytic model was ever fit to the SD profile of NGC 6397 extending to large projected radii⁹. Trager et al. (1995) estimated an effective (projected half-light) radius $R_{\mathrm{e}}=2\stackrel{\prime }{.}90$ for NGC 6397 using ground-based data and analyzing the SD profile in different apertures. Using Gaia DR2 data, de Boer et al. (2019) recently inferred $R_{\mathrm{e}}=3\stackrel{\prime }{.}35$ (which we infer from their 3D half-luminosity radius r_h, coupled with their assumed distance and with the ratio R_e/r_h ≃ 0.732 kindly estimated for us by M. Gieles). But their analysis was based on a fit to two dynamical models that both assumed a cored inner density profile, which clearly does not represent the inner SD profile of NGC 6397.

We thus had to estimate the SD profile ourselves. The cuspy profile led us to fit a Sérsic model (Sérsic 1963; Sersic 1968) to the distribution of projected radii since the Sérsic model has a shape parameter (Sérsic index) that makes it flexible to handle both cuspy and cored density profiles (e.g., Fig. A.1 of Vitral & Mamon 2020). Furthermore, there are simple analytical forms for the precise (but not exact) deprojection of the Sérsic model. In Appendix A, we extend the analysis of Vitral & Mamon (2020) of the precision of the different analytical approximations as a function of R_e and n, to even smaller values of R/R_e representative of our data. MAMPOSST-PM will then use the best-fit values of log R_e and n and their uncertainties as the respective mean and standard deviation of the Gaussian priors.

We required HST data to fit correctly the inner cusp, and we also needed Gaia data to fit correctly the outer envelope. We therefore fitted the GC surface density profile Σ_GC(R) together with the field star surface density Σ_FS (assumed constant) from the Gaia DR2 photometric data, supplemented by the HST data to better probe the very inner regions. The two datasets are not cleaned in any way, except for cuts in magnitude and maximum (minimum) projected radius. The HST data may not be radially complete (we did not have a parent sample with all stars including those without PM measurements), but this is the best we can do. Our approach does not remove the Milky Way field stars, but instead incorporates them in our SD fit. It allows us to keep Gaia stars with no PM information, which would be filtered out in a PM cut, even though their contribution to the surface density shape is relevant.

7.1.3. Stitching HST and Gaia

We stitched the HST and Gaia datasets as follows:

1. We cut the Gaia stars (GC stars plus interlopers) at G = 17, where we expect negligible mass segregation. We explored different magnitude limits, but moderate changes affected little our fits. However, using very faint magnitude limits (G = 19 or 20.5), we noticed a lack of continuity in the transition from the HST to the Gaia radial range, which appears to be a consequence of overcrowding in the cluster’s center (Arenou et al. 2018), where fainter magnitude stars are more difficultly detected.

2. We counted the number 𝒩_G of Gaia stars inside the annulus $2\stackrel{″}{.}7<R<100″$, corresponding to the region where our HST data seems complete (see Fig. 8).

3. We removed the Gaia stars in that annulus and added the 𝒩_G HST brightest stars (a tiny fraction of which may be interlopers) in that annulus, so we could correctly probe the inner surface density profile.

7.1.4. Surface density results

We performed maximum likelihood estimation (MLE) to fit log R_e and n to the distribution of projected radii, using a model SD profile with a Sérsic GC plus a uniform (constant SD) field star contribution, considering the data up to a maximum projected radius, R_max¹⁰. This maximum radius must be chosen large enough to capture the behavior of the SD profile in the outer envelope of the GC as well as the contribution of the field stars, but not too large to be overwhelmed by the field stars in the MLE fit.

We performed such MLE fits of the SD profile for 500 values of log-spaced R_max between $0\stackrel{°}{.}2$ and 1°. Figure 10 (left) shows a plateau for R_max in the range ∼1300″ − 1800″, where R_e, n and the field star SD are each indeed constant. At lower (resp. higher) R_max the field star SD fluctuates at lower (higher) values with corresponding higher (lower) and fluctuating R_e and n.

For better precision, we then performed 30 MCMC runs for log-spaced R_max between 1300 and 1800 arcsec and assigned the mean of best-fitted values and their uncertainties (green crosses and error bars of Fig. 10, respectively) as the global Sérsic radius and index means and uncertainties. The results are displayed in Fig. 10 (bottom). We found $R_{\mathrm{e}}=4\stackrel{\prime }{.}51\pm 0\stackrel{\prime }{.}36$ and n = 3.26 ± 0.23. Our effective radius is significantly higher than the estimates of Trager et al. (1995) and de Boer et al. (2019), which both fail to properly capture the cusp. Indeed, the Sérsic SD profile with the same index (n = 3.26) but with the Trager et al. (1995) effective radius (green curve in the right panel of Fig. 10) fails to reproduce the steep inner SD profile.

7.2. Multiple populations

Mass-orbit modeling codes such as MAMPOSST-PM assume a given population with a given SD profile, as well as given kinematics (i.e., p(v|R) for MAMPOSST-PM). Systems with multiple populations, each with their SD profiles and kinematics, should be analyzed jointly. MAMPOSST-PM can handle such multiple populations. In practice, we predetermine priors for the SD profile of each population, and then run MAMPOSST-PM with several tracer components.

There are two mechanisms that bring multiple populations in GCs. First, the formation of a GC may occur in several episodes, with different populations with their own chemistry and kinematics. Second, the short two-body relaxation time of GCs in general and especially of core-collapsed GCs such as NGC 6397, allows for energy exchanges between stars that are sufficient to drive them toward energy equipartition. This leads to mass segregation, where massive stars are confined to smaller radii and have lower typical velocities at a given radius. We examine these two mechanisms in turn.

7.2.1. Multiple chemical populations

Many GCs show signs of different chemical populations (e.g., Carretta et al. 2009a), in particular NGC 6397 (Carretta et al. 2009b; Milone et al. 2012b). Unfortunately, the separation of stars into different chemical populations in NGC 6397 was only performed with HST data limited to small projected radii. Such a chemical separation of stars is impossible with the Gaia data, for lack of sufficient wavebands to provide a large enough dimensionality of colors to distinguish chemical populations. Moreover, this requires cleaning the data for differential reddening. Therefore, the splitting of the stars of NGC 6397 into chemical populations, while possible, is beyond the scope of the present article.

Besides, Cordoni et al. (2020), who analyzed the spatial distributions and kinematics of seven GCs (other than NGC 6397) split by their two detected chemical populations, found that only the two GCs with the highest two-body relaxation times showed signs of different kinematics between their two chemical populations. If the two populations have roughly the same mass functions, then two-body relaxation should wash out any differences in their positions in projected phase space (Vesperini et al. 2013), except in their outer regions where two-body relaxation is incomplete. The very low relaxation time of NGC 6397 (600 Myr at its half-mass radius according to Baumgardt et al. 2019) thus suggests that its chemical populations should have similar kinematics, at least within its effective radius of $4\stackrel{\prime }{.}5$ as we found above.

7.2.2. Mass segregation

We explored the range of masses of the stars in our datasets by comparing synthetic CMDs to our observed ones. Following Carballo-Bello et al. (2012), we generated a CMD with a set of isochrones from the PARSEC (Marigo et al. 2017; Bressan et al. 2012; Pastorelli et al. 2019 and references therein) software. We provided n_inTPC = 15 since Marigo et al. (2017) mention this suffices to recover the main details of thermal pulse cycle variations in the evolutionary tracks. We also assigned η_Reimers = 0.477, which is the median value of the Reimers coefficient (Reimers 1975) among clusters, provided by McDonald & Zijlstra (2015); in their Fig. 4, they showed that NGC 6397 presents a value of η_Reimers corresponding to this value within ∼1σ. We adopted an age of 12.87 Gyr and a metallicity of [Fe/H] = −1.54 from Marín-Franch et al. (2009), while adopting our distance of 2.39 kpc (Sect. 3.3) and a Galactic extinction A_V = 0.50 obtained from E(B − V) = 0.160 (Schlafly & Finkbeiner 2011)¹¹ and R_V = 3.1. We note that this value of extinction is 14% lower than E(B − V) = 0.188 found in the older reddening maps of Schlegel et al. (1998), but consistent with A_V = 0.52 ± 0.02 found by Valcin et al. (2020) despite their larger distance of 2.67 pc (Table 1).

The left panel of Fig. 11 shows that the PARSEC model fits well the CMD of the cleaned subsample of HST stars. The mass ratio among the heaviest and the lightest stars in our cleaned HST sample is 2.57, sufficiently high to worry about the possibility of mass segregation in our cleaned HST sample.

In fact, mass segregation is visible in NGC 6397 (Heyl et al. 2012; Martinazzi et al. 2014). Analyzing HST stars in the range $3\prime <R<7\stackrel{\prime }{.}5$, Heyl et al. (2012) found mass segregation of main-sequence stars: brighter stars show 4% lower median projected radius and 12% lower median PM moduli than fainter stars. In a subsequent study by the same team, Goldsbury et al. (2013) showed that two characteristic radii scale as a M^γ, with γ = −1.0 ± 0.1 for NGC 6397. Martinazzi et al. (2014) found that the mean mass of main-sequence stars drops with physical radius (after deprojection) from > 0.7 M_⊙ for r <  10″ to ≃0.56 M_⊙ for r >  100″ (after correcting for the completeness with magnitude, estimating it by adding artificial PSF-convolved stars to the images), thus a 20% effect.

MAMPOSST-PM is able to treat multiple stellar populations together. We therefore also performed SD fits with two populations of stars, as we explain below.

The magnitude threshold we chose to separate the two population of stars was based on the analysis of Heyl et al. (2012): Their Figs. 7 and 12 indicate that NGC 6397 main-sequence stars present a (small) variation of radial distribution and velocity dispersion profiles, respectively, at a magnitude of F814W ∼ 18.75, which corresponds to F606W = 19.76. We thus used this limit to divide our cleaned subset of 8255 stars into two populations, one with brighter magnitudes (6527 heavier stars) and the other with fainter ones (1728 less massive stars).

Unfortunately, we could not perform a surface density fit as robust as for the single population case because (1) our data was considerably incomplete at higher magnitudes (the fainter subset) and (2) as mentioned in Sect. 7.1.3, the HST plus Gaia stacked subset presented discontinuous trends when allowing stars with fainter magnitudes (G >  17). We therefore let MAMPOSST-PM fit the SD profile of each population from the kinematics only, more precisely from the conditional probabilities p(v|R). We adopted Gaussian priors for log R_e, with mean equal to that found by our previous SD fit of the single population (Sect. 7.1) and a wide (0.2 dex) uncertainty. We also considered Gaussian priors for the Sérsic indices, but with lower uncertainties to avoid a degeneracy in its marginal distribution.

We helped MAMPOSST-PM by providing narrow Gaussian mass priors for each population. We derived mass fractions for each population using the power-law main-sequence stellar mass function of slope α = −0.52¹² together with the mass limits of our subsets. Indeed, a power-law relation dN/dm ∝ m^α implies a total stellar mass in the range of stellar masses (m₁, m₂) of

$$\begin{array}{c}\hfill M_{\mathrm{total}}\propto {\displaystyle {\int }_{m_1}^{m_2}m\frac{\mathrm{d}N}{\mathrm{d}m}\mathrm{d}m=\frac{m_2^{2+\alpha }-m_1^{2+\alpha }}{2+\alpha },}\end{array}$$(22)

and thus derive

$$\begin{array}{c}\hfill \frac{M_{\mathrm{bright}}}{M_{\mathrm{faint}}}=\frac{m_{\mathrm{bright}}^{2+\alpha }-m_{\mathrm{cut}}^{2+\alpha }}{m_{\mathrm{cut}}^{2+\alpha }-m_{\mathrm{faint}}^{2+\alpha }},\end{array}$$(23)

where M_bright and M_faint are the masses of the brighter and fainter populations, respectively and m_bright, m_faint and m_cut are the respective highest, lowest and two-population threshold masses of the global subset. With m_bright = 0.77 M_⊙, m_cut = 0.51 M_⊙, and m_faint = 0.25 M_⊙ (right panel of Fig. 11), Eq. (23) yields a bright mass fraction of 0.56. Since our main mass estimates of NGC 6397’s mass with one single population were centered around 10⁵ M_⊙, we passed logarithmic Gaussian priors to each population, centered at log(M/M_⊙) = 4.7, with standard deviation of 0.05.

7.2.3. MAMPOSST-PM assumptions and priors

To summarize, MAMPOSST-PM assumes that NGC 6397, lying at a distance of 2.39 kpc, with its center at the position found by Goldsbury et al. (2010), is a spherical self-gravitating system with no internal streaming motions (rotation, expansion/collapse, etc.). It also assumes that the local velocity distribution function is an ellipsoidal 3D Gaussian aligned with the spherical coordinates.

By relying on the Jeans equation, MAMPOSST-PM assumes that stars are test particles orbiting the gravitational potential, and that they do not interact with one another. While binary stars have their own motions and violate the test particle hypothesis, these are in great part filtered out of our data. However, the low two-body relaxation time of NGC 6397 means that star-star interactions are not negligible, which violates the test particle assumption.

We now summarize our adopted priors. The surface density profile of GC stars is assumed to follow the Sérsic law. We adopted Gaussian priors for log R_e and for n (see Sect. 7.1) for both single- and two-population analyses.

Following Sect. 5.2.3, we adopted Gaussian priors on the bulk LOS velocity and PM of the GC and Pearson VII priors on the modulus of the bulk PM of the field stars. The field star LOS distribution is Gaussian with a wide (1.7 dex) uncertainty. The log ratio of field stars to GC stars is flat between −5 and −2.

Our standard runs assume isotropic velocities. We also ran cases with anisotropic velocities, and following Mamon et al. (2019) we used flat priors on β_sym = β/(1 − β/2), which varies from −2 for circular orbits to 0 for isotropic orbits and +2 for radial orbits, and with β_sym → β for |β| ≪ 1. Our flat priors were limited to −1.9 <  β_sym <  1.9, where the two extremes correspond to β = −38 and β = 0.97, respectively (see Eq. (8)).

We also adopted flat priors on the GC stellar log mass (3 to 6) for the single-population runs and Gaussian priors with equal means and wide dispersions for the two-population runs (Sect. 7.2.2). The (sub-)cluster of unseen objects is modeled with different density models with a wide Gaussian (1 dex) prior on the scale radius uncertainty, assuming isotropic velocities.

8. Analysis

8.1. Marginal distributions and covariances

We explored parameter space to determine marginal distributions and parameter covariances using the same Markov chain Monte Carlo (MCMC) approach as Mamon et al. (2013, 2019), see Sect. 4.7 of Mamon et al. (2019). This makes use of the public Fortran COSMOMC code (Lewis & Bridle 2002)¹³. In particular, we used 6 MCMC chains run in parallel and stopped the exploration of parameter space after one of the chains reached a number of steps N_steps = 10 000 N_free, where N_free is the number of free parameters of the model. The only difference is that we now discard the 3000 N_free first steps of each MCMC chain (burn-in phase) instead of the first 2000 N_free as Mamon et al. (2019) previously did.

8.2. Bayesian information

There are several ways to compare the different results of MAMPOSST-PM using different priors. The simplest would be to compare the log-likelihoods. But there is then the risk of over-fitting (under-fitting) the data when using too many (few) free parameters. We considered two model selection criteria to distinguish our different models and priors. We first considered the corrected Akaike Information Criterion (Sugiura 1978)

$$\begin{array}{c}\hfill \mathrm{AICc}=\mathrm{AIC}+2\frac{N_{\mathrm{free}}(1+N_{\mathrm{free}})}{N_{\mathrm{data}}-N_{\mathrm{free}}-1},\end{array}$$(24)

where AIC is the original Akaike Information Criterion (Akaike 1973)

$$\begin{array}{c}\hfill \mathrm{AIC}=-2ln{\mathcal{L}}_{\mathrm{MLE}}+2N_{\mathrm{free}},\end{array}$$(25)

and where ℒ_MLE is highest likelihood found when exploring the parameter space, N_free is the number of free parameters, and N_data the number of data points¹⁴. We also considered the Bayes Information Criterion (BIC, Schwarz 1978):

$$\begin{array}{c}\hfill \mathrm{BIC}=-2ln{\mathcal{L}}_{\mathrm{MLE}}+N_{\mathrm{free}}ln{N}_{\mathrm{data}}.\end{array}$$(26)

Given that the relative likelihood (given the data) of one model relative to a reference one is

$$\begin{array}{c}\hfill exp(-\frac{\mathrm{AIC}-{\mathrm{AIC}}_{\mathrm{ref}}}{2})\end{array}$$(27)

(Akaike 1983), a model with a higher AICc (or AIC) relative to a reference model can be ruled out with 95% confidence if AICc >  AICc_ref + 6. BIC follows the analogous to Eq. (27) for the relative likelihood (Kass & Rafferty 1995):

$$\begin{array}{c}\hfill exp(-\frac{\mathrm{BIC}-{\mathrm{BIC}}_{\mathrm{ref}}}{2}),\end{array}$$(28)

so a model with a higher BIC relative to a reference model can be ruled out with 95% confidence if BIC >  BIC_ref + 6.

Since BIC penalizes extra parameters more (factor lnN_data ≃ 9.0) than AICc (i.e., by a factor two for our dataset), BIC effectively prefers simpler models than does AICc. BIC is well suited for situations where the true model is among the tested ones, while AIC(c) is more robust in the opposite case (Burnham & Anderson 2002). In practice, we do not expect any of our models to be true: for example, there is no expectation that the surface density of stars in a GC should precisely follow a Sérsic model. To paraphrase the statistician George Box: “All models are wrong, but some are useful” (Box 1979)¹⁵.

We therefore present both AICc and BIC in our results. However, we give preference to AICc whenever having to decide between different models.

9. Mass-orbit modeling results

9.1. MAMPOSST-PM runs

Table 3 displays our standard runs, along with some robustness tests on velocity isotropy, mass segregation and presence of central components. Our runs have come with 4 variations for the inner mass distribution: (1) no IMBH nor (sub)-cluster of unresolved objects (CUO, see Sect. 2.2), (2) IMBH with no CUO, (3) CUO without an IMBH, and (4) both IMBH and CUO.

The left panel of Fig. 12 shows the marginal distributions and covariances from the MCMC output of model 1 (IMBH without CUO and with isotropic stars), for a selection of missing parameters (missing are the bulk GC motion and the log ratio of field stars to GC members). The MCMC ran well in the sense that the MLE values (red arrows) are close to the peaks of the marginal distributions (blue shaded histograms). One notices that higher IMBH masses are obtained with higher GC scale radii, independently of the Sérsic index. The reason is that, at given total stellar GC mass, more concentrated stars will lead to more mass in the inner regions, meaning less IMBH mass. This is the opposite of what van der Marel & Anderson (2010) proposed, but consistent with what Baumgardt (2017) found for NGC 6397.

The right panel of Fig. 12 shows the same for model 10 (CUO without IMBH and with isotropic stars). Again, the MLE values are coincident with the peaks of the marginal distributions. But surprisingly, the MCMC favors a GC scale radius at the upper edge of the 3σ prior obtained from the SD profile fit of Sect. 7.1.4. We will discuss this below. The mass of the CUO increases with the GC and CUO scale radii, and decreases with the GC Sérsic index (because low index leads to flatter SD profiles, which are somewhat analog to wider SD profiles). Interestingly, the scale radius of the CUO population is small (7″), as we will discuss in Sect. 9.4.1.

9.2. Velocity anisotropy

The combination of HST data properly probing the inner regions of NGC 6397 with the Gaia data probing the outer regions allows us to estimate the velocity anisotropy across the cluster. We ran MAMPOSST-PM using different priors on the anisotropy. Our standard prior has isotropic velocities throughout the GC. Our other priors assume the gOM anisotropy model (we found that the softer varying Tiret et al. 2007 model performs almost as well, but not better).

The highest likelihoods (lowest −Δlnℒ_max) between the first five models are achieved by model 2 (everything free), but it did not converge (R⁻¹ = 0.07 >  0.02). In particular, model 2 leads to quasi-isotropic inner and outer velocities, with a narrow constraint for β₀ and a quite wide one for β_∞. The transition radius, r_β, is very poorly constrained. Compared to model 2, the isotropic model 1 is very strongly preferred by BIC (according to Eq. (28) for ΔBIC = 26) and marginally so (92% confidence according to Eq. (27) for ΔAIC = 4.87) by AICc. Model 3 with the anisotropy radius tied to the effective radius of the SD profile (TAND) is strongly preferred by BIC Bayesian evidence (with respect to model 2), but not with AICc. Interestingly, model 3 leads to negligible anisotropy in the center with narrow uncertainty. But it points to mildly tangential anisotropy at r = r_out = 8′, which corresponds to our maximum projected radius R_max. Forcing central isotropy with TAND (model 4) leads to less tangential outer velocities. Conversely, forcing outer isotropy with TAND (model 5) leads to strongly well constrained isotropic inner velocities.

We also ran a CUO case, where the stars had free inner and outer anisotropy, but with TAND. Compared to its isotropic analog (model 10), model 13 with a single population has a slightly worse likelihood (but did not fully converge), and the inner and outer anisotropies are consistent with isotropic orbits (β₀ = −0.02 ± 0.05, β_out = −0.03 ± 0.05). An analogous run for the two-population CUO model also yields close to isotropic orbits for the bright population (β₀ = 0.03 ± 0.05, β_out = 0.00 ± 0.12). The faint population has inner velocities consistent with isotropy (β₀ = −0.1 ± 0.1). On the other hand, its outer anisotropy is poorly constrained (${\beta }_{\mathrm{out}}=0.1_{-0.7}^{+0.2}$) because of the near complete lack of faint Gaia stars with decent PM errors, hence we are limited to the inner regions. Still, the outer velocities are again consistent with isotropy.

In summary, there is strong evidence for isotropy with BIC, but weaker evidence with AICc. However, even the anisotropic runs produce anisotropy profiles that are very close to isotropic throughout. We thus conclude that the visible stars in NGC 6397 have quasi-isotropic orbits, at least for the stars brighter than F606W = 19.76 and up to R = 8′. We therefore adopted isotropic orbits as our standard when investigating other quantities.

9.3. Intermediate mass black hole

We tested the scenario with an IMBH and no CUO (models 1 to 6). Among models 1 to 6, the most likely one and very strongly favored by AICc (model 6, which is the analog of model 1, but with 2 populations for the GC stars) yields an IMBH mass of ${511}_{-207}^{+158}{M}_{\odot }$, while model 1, which is weakly favored by BIC over model 6, yields an IMBH mass of ${658}_{-338}^{+70}{M}_{\odot }$. Both models 1 and 6 indicate an IMBH mass above 200 M_⊙ at 95% confidence. Furthermore, AICc (resp. BIC) indicates very strong (resp. quite strong) evidence for the presence of an IMBH in the absence of a central diffuse component (i.e., comparing IMBH models 1 and 6 to models 7 and 8, respectively).

However, both AICc and BIC indicate strong evidence against the IMBH hypothesis in comparison with the presence of a CUO, with differences of 15.9 in AICc and 8.9 in BIC between isotropic, single-population models 1 and 10. With two-population GC stars, we have differences of 12.2 in AICc between isotropic models 6 and 14. But with BIC, the difference is only 5.2, leading to only moderately strong evidence (92% confidence) of favoring the CUO over the IMBH. But given our preference for AICc for the complex physics of GC kinematics with imperfect models for the surface density profile for example (see Sect. 8), we conclude that the evidence is strong in favor of a dark component that is diffuse for the two-population model. Hence, the unseen inner matter of NGC 6397 is very likely to be diffuse. We therefore now investigate the CUO model in more detail.

9.4. Inner subcluster of unresolved objects (CUO)

9.4.1. CUO density profile

If the dark component is diffuse as a CUO instead of a singular IMBH, we first need to measure its extent. MAMPOSST-PM can provide constraints on the shape and scale radius of the density profile of the unseen population, through the sole use of the conditional probabilities of velocity at given projected radius, p(v|R), without directly fitting the distribution of projected radii.

We first assumed a Plummer (1911) model for the CUO with the same effective radius as that of the GC stars ($4\stackrel{\prime }{.}5$), but with a wide standard deviation (1 dex) for the Gaussian prior on log scale radius. Interestingly, as seen in Fig. 13, MAMPOSST-PM converged to r₋₂ = 7″, thus an effective (half-projected number) radius of only $4\stackrel{″}{.}3$ for the CUO¹⁶, which is 60 times lower than the median of the prior, thus confirming our suspicion that the CUO might indeed be significantly more concentrated than the main-sequence stars.

We then ran MAMPOSST-PM for three models with a much smaller scale radius prior for the CUO (30″), to allow for more accurate fits. We used three different CUO density models: a Sérsic model, a Plummer (1911) model, and a Hernquist (1990) model. All three led to very small scale radii, which convert to effective (half projected number) radii $R_{\mathrm{e}}=9\stackrel{″}{.}6$, $17\stackrel{″}{.}4$ and $4\stackrel{″}{.}4$ for the Sérsic, Hernquist and Plummer models, respectively¹⁷.

The Sérsic model (12) is the most likely one, closely followed by the Plummer model (10). AICc has a weak preference for the Plummer model, while BIC prefers the Plummer model, with strong evidence against the more complex Sérsic model, but weak evidence against the Hernquist model (11). Model 12 with the Sérsic CUO density profile produced a low Sérsic index of n = 0.92, which leads to a shallow inner slope that is not too different from the zero slope of the inner Plummer density profile (see Fig. A.1 from Vitral & Mamon 2020). The steeper inner profile of the Hernquist density profile makes it less similar to the n = 0.92 Sérsic model than is the Plummer profile, as confirmed by the Hernquist model 11 producing the lowest likelihood (highest −lnℒ) among the three models.

In summary, it is hard to distinguish which is the best density model for the CUO scenario. There is weak evidence for a shallow slope. We adopted the Plummer model given that it is the preferred of the three density models for both AICc and BIC, albeit with weak evidence for both.

9.4.2. Presence of an IMBH in addition to the CUO

One may ask whether the center of NGC 6397 can host both an IMBH and a diffuse dark component. Model 9 contains both, but it is somewhat less likely than model 10, and in comparison it is strongly disfavored by BIC, although only weakly disfavored by AICc. We note that the Plummer model used for the CUO is the one that best distinguishes the CUO from a possible additional IMBH. Moreover, the recovered mass of the additional IMBH is so small ${42}_{-26}^{+92}{M}_{\odot }$ that it can no longer be called an IMBH.

9.5. Two-mass populations

Table 3 also displays the MAMPOSST-PM results of 4 models with two populations, split by apparent magnitude, hence by mass (Sect. 7.2.2). Model 14, with CUO but no IMBH, has the highest likelihood of the three models considering isotropy. It is very strongly favored with AICc over models 6 (with IMBH but no CUO) and 8 (no IMBH nor CUO). It is also preferred by BIC, which strongly favors it over model 8, but only marginally favors it over model 6 (ΔBIC = 5.16, leading to 92% confidence in preferring model 14 over model 6, according to Eq. (28)). This preference of the CUO model over the IMBH model resembles that found for the single population (Sect. 9.3).

But there are differences in the MAMPOSST-PM results between single-population and two-population models, both in their Bayesian evidence and in their best-fit parameters. Indeed, the two-population model 14 shows the highest likelihood and AICc evidence of the first fourteen models listed in Table 3. In particular, comparing two-population models to their single-population equivalents, there is strong AICc evidence (ΔAICc = 13.7) favoring two-population model 14 than 1-population model 10. But there is strong BIC evidence (ΔBIC = 7.3) the other way, with model 10 displaying the best BIC evidence of all fourteen models listed in Table 3.

The reader may note that the GC stellar masses summed over the one or two populations are 9% greater in the two-population runs for CUO without IMBH scenarios (model 14 versus model 10) and 13% higher for IMBH without CUO scenarios (model 6 versus model 1). Moreover, the IMBH mass is 22% lower in the two-population model 6 relative to the single-population model 1. Given that the standard deviations in the log GC mass (single or sum of Bright and Faint) are less than 0.03 dex, the difference in the means of log M_GC of 0.086 dex is highly significant for samples of order of 10⁵ points, as we checked with a Student t test. Similarly, the standard deviations in log M_IMBH are of order 0.2 dex; thus the difference of 0.037 dex in their means is again highly significant, as we also checked with the Student t test.

MAMPOSST-PM yields interesting results on the differences between the bright and faint populations. First, the two-population runs can be tested for the respective masses in each. Despite our prior of equal masses for each, MAMPOSST-PM returns best-fit bright population mass fractions of 0.49, 0.44, and 0.57 for models 6 (IMBH without CUO), 8 (no IMBH nor CUO), and 14 (CUO without IMBH). Only model 14 has a bright fraction close to the expected value of 0.56 (Sect. 7.2.2).

Secondly, in all three two-population models, the brighter population has a much lower scale radius than its fainter counterpart, by factors of 2.5, 2.6 and 2.1 for models 6, 8, and 14, respectively. These lower scale radii for the bright population are highly statistically significant. Indeed, the fractions of MCMC chain elements leading to higher scale radius of the brighter population are less than 0.04% for models 6, 8, and 14. Therefore, MAMPOSST-PM is able to find very strong kinematic signatures of luminosity (hence mass) segregation, by fitting p(v|R) with the same priors on the scale radii of the two populations, without directly fitting the distribution of projected radii.

Recall that we used the same Gaussian priors on scale radius for both populations. As can be seen in Fig. 14, these priors may have been too narrow, and MAMPOSST-PM may have thus under-estimated the differences in the scale radii of the bright versus faint populations.

Figure 15 illustrates the quality of the 3 classes of models in reproducing the observed velocity dispersion profiles. Model 8, with no additional dark component, clearly underestimates the velocity dispersions below 10″. Model 6, with an IMBH, overestimates it below 4″. Finally, model 14, with a CUO, does best (the non-perfect match appears to be caused by sampling noise).

10. Conclusions and discussion

10.1. Main results

We ran the MAMPOSST-PM code on a combined set of PMs from HST and Gaia, as well as with LOS velocities from MUSE. Our results indicate that the GC star velocities are close to isotropic out to 2 effective radii. Our models with inner central components (IMBH or diffuse subcluster of unresolved objects, CUO) are strongly preferred over models without any. Models where the central mass is in the form of an IMBH find favored masses of 500 − 650 M_⊙, with a 5th lower percentile of 200 M_⊙. But we found better likelihoods and Bayesian evidence for a diffuse central unseen mass component (CUO), with effective (projected) radius of order of 2.5 to 5″ and mass of order 1000 to 2000 M_⊙.

10.2. Robustness

We tested the robustness of our results for several variations in our assumptions and our cuts to the data. These are: Single versus two stellar populations, the adopted SD profile, the minimum allowed PM error, and the HST quality flags.

First, the preference for the CUO component is relatively robust to our analysis with a single bright population or two populations. It is also quite robust to our assumptions on velocity anisotropy from free inner and outer values to isotropic. The differences in AICc and BIC of IMBH relative to CUO are +16 and +9 for single-population isotropic, +15 and +8 for single-population, free inner and outer anisotropy but fixing the anisotropy transition radius to the density scale radius (TAND), and +12 and +5 for double-population isotropic. In other words, CUO is very strongly preferred to IMBH by both AICc and BIC. We also ran an isotropic model restricting the sample to stars with F606W and G magnitudes brighter than 17.5, to limit the effects of possible mass segregation of these bright stars with the fainter, lower-mass stars. We found, again, that the CUO scenario is strongly preferred over the IMBH one (ΔAICc = 7.0).

The preference for CUO versus IMBH is also robust to the choice of the density profile. Indeed, adopting a Sérsic profile, but fixing the effective radius to $2\stackrel{\prime }{.}9$ as found by Trager et al. (1995), produces a worse match to the kinematic data. However, it strongly prefers the CUO model to the IMBH one, but provides much weaker constraints on the central masses. Similarly, forcing an inner cored profile to the luminous mass, with a Plummer model, keeping the same effective radius, the kinematics indicate much worse AICc compared to Sérsic. The AICc then strongly prefers anisotropic models, for which it still strongly prefers the CUO option to the IMBH one.

In addition, the preference for a diffuse central mass (CUO) versus an IMBH is also robust to the choice of maximum allowed PM error. We obtain the same ΔAICc = 12 for both ϵ_μ <  0.276 mas yr⁻¹ and our standard choice of ϵ_μ <  0.197 mas yr⁻¹, as well as the same ΔBIC = 5. On the other hand, to ϵ_μ <  0.141 mas yr⁻¹, yields ΔAICc = 6.5 instead, presumably because of the small dataset satisfying this low maximum allowed PM error.

Finally, we also tested the robustness of our conclusions when also applying HST PM quality flags, such as the reduced χ² flag. Limiting our sample to stars satisfying $({\chi }_{{\mu }_{\alpha }\text{\hspace{0.05em} cos}\delta }^2+{\chi }_{{\mu }_{\delta }}^2)/2<1.25$, which is the threshold chosen by Bellini et al. (2014) for NGC 7078, we noticed that the fitted parameters agreed within 1σ to the ones in Table 3. Furthermore, the CUO scenario was again favored over the IMBH one with moderately strong AICc evidence (ΔAICc = 4.91, i.e., 91% confidence) instead of the strong AICc evidence observed for the subset without this χ² filter. Moreover, Bellini et al. (2014) argue that no quality parameter filters are required for nearby GCs with multiple observations, such as NGC 6397, when the PM dispersion depends little on them, as we indeed checked for both the reduced χ² and for their QFIT parameter. Our conclusions are thus robust with respect to the HST data quality flags.

10.3. Effects of the dataset

It is useful to check how efficient are our different datasets in reaching our conclusions. Such a comparison will be useful for future mass-orbit modeling of other GCs, for example in the absence of HST and LOS data.

Table 4 displays a summary of the masses and CUO size from single-component isotropic runs using different datasets. We used the same SD profile for all samples.

First, with MUSE data alone, we cannot constrain the mass of the IMBH. In contrast, Kamann et al. (2016) found M = 600 ± 200 M_⊙ for the IMBH. The difference can be explained by the much more liberal cut in velocity errors used by Kamann et al. (2016), 5 km s⁻¹, compared to our cut at 2.2 km s⁻¹, which led them to have 4608 stars compared to our 528 stars. Our small MUSE sample prevents us from constraining the CUO mass and size.

This is also true using Gaia data alone. It is not surprising, given that we have only 152 Gaia stars within the inner 100″ and none within the inner 10″. This low number is caused by our cut on PM errors. Later Gaia releases will have lower PM errors, and will thus have larger numbers of inner stars to use for mass-orbit modeling.

Moreover, a subset with just HST and MUSE, extending up to $2\stackrel{\prime }{.}5$ is not able to constrain the outer isotropy indicated by Gaia DR2. Indeed, a run similar to model 13 with this dataset yielded $\beta (r=8\prime )=-1.0_{-0.3}^{+0.5}$, which is a significantly tangential anisotropy, in contrast with $\beta (8\prime )=-0.{03}_{-0.14}^{+0.11}$ when Gaia data is also considered.

In summary, Gaia offers much improvement to HST(+MUSE) data, by (1) allowing a better estimation of the surface density profile, hence of the 3D density profile, and (2) allow measuring velocity anisotropy well beyond the effective radius. With the prior knowledge of the surface density profile and velocity anisotropy profile, the addition of Gaia provides somewhat better constrained masses for the inner component (IMBH or CUO). On the other hand, Gaia data by itself, has much too few stars with sufficiently accurate PMs at low projected radii to probe, the IMBH or CUO mass, and even less the CUO size. This highlights the requirement of combining HST and Gaia for proper mass-orbit modeling of the PM of stars in GCs.

10.4. Velocity anisotropy

Our combination of HST data probing the inner regions of NGC 6397 and Gaia data probing the outer regions, allowed us to obtain a wide-range view of the variation of the orbital anisotropy with radius. Bayesian evidence favors isotropy throughout the cluster. This confirms the projected isotropy previously determined by Heyl et al. (2012) from $3\stackrel{\prime }{.}5$ to 7′ using HST (with 10% individual errors), as well as Watkins et al. (2015a): σ_POSt = 0.98 σ_POSr, up to R_e/10 using HST, and by Jindal et al. (2019) at larger radii up to 12′≃3 R_e using Gaia, with signs of radial motions at R >  8′.

But a constant projected isotropy can hide radial variations of the 3D anisotropy β(r). Figure 16 displays the first known constraints on the velocity anisotropy profile of NGC 6397, which we obtained from our MAMPOSST-PM analyses from 1″ to 8′. Indeed, it is highly consistent with isotropic orbits throughout, as seen in Fig. 16. Our parametric representation of β(r) prevents us to see a possible sharp up- or down-turn in this radial range. Nevertheless, MAMPOSST-PM, like other mass-orbit modeling algorithms, probes physical radii beyond the maximum projected radius. This suggests that the orbits of stars at two effective radii are not much disturbed by the Milky Way, despite the differences in the two components of the mean POS velocity beyond R = 8′.

Only two previous studies have measured velocity anisotropy profiles in GCs. Cluster ω Cen was modeled by van der Marel & Anderson (2010), who fit the PM dispersions measured in bins of projected radius. They found a gOM anisotropy profile with β₀ = 0.13 ± 0.02 and β_∞ = −0.52 ± 0.22, thus slightly radial in the center and slightly tangential in the outer regions. Messier 15 was analyzed by den Brok et al. (2014), who used an updated version of JAM (Cappellari 2008). They found a gOM anisotropy profile with β₀ = −0.21 ± 0.30 and β_∞ = 0.015 ± 0.12, hence consistent with isotropic at all radii.

In comparison, the present study provides the first Bayesian mass-orbit modeling of discrete GC data. Our model 13, as well as the bright population of model 15, both yield inner anisotropies consistent with zero with ≈ ± 0.05 uncertainty (Table 3). Both yield anisotropies at 8′ of |β|≤0.03 with uncertainties ≈0.14 (Table 3).

The velocity isotropy of NGC 6397 at fairly large radii appears to contradict models where core-collapsed GCs have isotropic inner motions and very radial outer motions (e.g., Takahashi 1995). But such models are for isolated GCs, whereas NGC 6397 has crossed the Galactic disk many times and as recently as 4 Myr ago (Sect. 2.2). GC stars are perturbed by tidal forces from the Milky Way (disk, bulge and halo, which sum up to a total gravitational potential not far from spherical, see Fig. 2.19 of Binney & Tremaine 2008). Indeed, we noticed a difference in the mean POS velocities between the radial and tangential components beyond 8′ (top panel of Fig. 6), suggesting that the outer parts of NGC 6397 are out of equilibrium.

While the tidal field on an embedded GC inside a galaxy should be compressive and radial (Dekel et al. 2003), it is now understood that this does not lead to radial orbits in the outer envelopes of GCs, where the stars are most susceptible to tidal perturbations. Indeed, comparing stars of same binding energy, those that are on radial orbits spend more time at larger radii and feel stronger tidal forces, while those that are on tangential orbits are less tidally perturbed (see Giersz & Heggie 1997). This reduces the radial outer anisotropy caused by two-body relaxation. This is confirmed by the isotropy at the half-mass radius of GCs with short half-mass relaxation times determined with mass-orbit modeling by Watkins et al. (2015a), who also found increasingly radial orbits at r_h for GCs with increasingly longer half-mass relaxation times.

This is also consistent with N-body simulations of the internal motions of a GC subject to the tidal field of a point-mass galaxy (Tiongco et al. 2016; Zocchi et al. 2016) or of a realistic galactic potential by itself or within an infalling dwarf galaxy (Bianchini et al. 2017). These simulations indicate that, after the establishment of radial outer anisotropy at core collapse, the outer velocity anisotropy becomes less radial thanks to the tidal field of the Milky Way, especially when the tidal field is strong, as measured by small values of the ratio of Jacobi to half-mass radii, r_J/r_h, also called the filling factor. One can compare these simulations in detail to NGC 6397. We estimate the filling factor for NGC 6397 using its estimated pericentric and apocentric radii of 2.9 and 6.6 kpc (Gaia Collaboration 2018a) and the Milky Way mass profile of Cautun et al. (2020), which produces r_J = 31 to 51 pc, hence r_h/r_J = 0.08 to 0.14, assuming R_e/r_h = 0.74¹⁸. The single simulation of Zocchi et al. (2016) reaches this range of r_h/r_J for snapshots 3 to 11, for which β(r_h) decreases from 0.13 to −0.02, as the tides preferentially remove the radial orbits, leading to β(r_h) = 0.00 to 0.08 for r_h/r_J = 0.12, decreasing to more tangential values for higher filling factors. Similarly, among the MW-only simulations of Bianchini et al. (2017), who also included a realistic initial mass function and stellar evolution, only one of them has a filling factor at 10 Gyr that matches the one deduced above for NGC 6397, and this simulation produced β(r_h) = 0.075. So both studies would predict that a core-collapsed cluster like NGC 6397, subject to the tidal field of the Milky Way, should recover quasi-isotropic orbits at r_h, consistent with what we found. This needs to be checked with more refined N-body simulations that consider elongated orbits of a GC around its galaxy and possibly the incorporation of binaries at the start of the simulation.

There are several additional possible mechanisms for producing outer isotropic orbits. First, the passage through the differentially rotating Galactic disk could lead to exchanges of angular momentum leading to tangential orbits, in particular for the outer less bound stars. But on the scale of a GC, the differences in disk rotation velocity must be minute. Second, the stars in the outer parts of the GC must suffer violent relaxation during passages through the Galactic disk, but the passages are so short (a few Myr) that there may not be enough time for violent relaxation to effectively perturb the GC stars. Third, solenoidal modes of turbulence on Galactic gas clouds generate angular momentum within these clouds and within the stars formed therein, which in turn can transfer this angular momentum to the passing GC. But again, the short duration of these passages limits the amount of angular momentum that the GC stars can acquire. Therefore, the outer isotropy appears to be produced by the tidal field of the Milky Way preferentially removing those stars that were on radial orbits.

10.5. Intermediate mass black hole

Our IMBH mass estimates between 500 and 600 M_⊙, with uncertainties of ≈200 M_⊙, are consistent with those of Kamann et al. (2016) (M_IMBH = 600 ± 200 M_⊙) and Tremou et al. (2018) (M_IMBH <  610 M_⊙). Our 95% confidence lower limits (i.e., 5th percentiles) are 201 M_⊙ for both models 1 (single population) and 6 (two-population).

But if IMBHs grow by BH mergers, one expects that the gravitational radiation emitted during these mergers are anisotropic and lead to substantial recoil of the BH remnant of the merger (Peres 1962). Such recoils are strong enough for 70% of IMBHs of mass < 1000 M_⊙ to escape a GC of escape velocity 50 km s⁻¹ over time (Holley-Bockelmann et al. 2008). We can estimate the escape velocity from the GC center, $v_{\mathrm{esc}}=\sqrt{-2\mathrm{\Phi }(0)}$, for our best IMBH model 6, given that the central potential is Φ₀ = −(2/π) bⁿ Γ(n)/Γ(2n) G M/R_e for the Sérsic model (Ciotti 1991), where b(n) follows the relation given in Ciotti & Bertin (1999). Model 6 then has v_esc = 20.9 km s⁻¹ and 13.1 km s⁻¹ for the Bright and Faint components, respectively. Summing these two terms quadratically leads to a global escape velocity of v_esc = 24.7 km s⁻¹ for model 6. This escape velocity is lower than the 50 km s⁻¹ assumed by Holley-Bockelmann et al. (2008), which should lead to an even much higher fraction of escaping IMBHs of mass 1000 M_⊙. And our 95% confidence upper limits on the IMBH mass (without a CUO) are lower, again implying that a yet higher fraction of IMBHs escape their GC. A simpler model by Holley-Bockelmann et al. (2008) indicates that 25 IMBH+BH mergers should lead to 85% of 500 M_⊙ IMBHs escaping their GC of escape velocity 17 km s⁻¹. This weakens our confidence in the central IMBH scenario for NGC 6397.

10.6. Cluster of unresolved objects

The scale radius of the CUO in NGC 6397 is consistent with a feature in the SD profile of NGC 6397. Indeed, as seen in the right panel of Fig. 10, the SD profile suggests a separation between two populations at R = 10″, consistent with the CUO effective radius of a $2\stackrel{″}{.}5$ to 5″.

We now discuss the nature of the subcluster of unresolved objects that we found in the center of NGC 6397. The much smaller scale radius for the CUO indicates that the objects must be more massive than the stars that we studied, thus more massive than m_bright = 0.77 M_⊙ (Sect. 7.2.2). Indeed, such unresolved massive objects would sink to the center by dynamical friction Chandrasekhar (1943). The core-collapsed state of NGC 6397 is a strong argument in favor of a population of heavier stars in the cluster’s inner regions, as it suggests a short two-body relaxation time, hence mass segregation. These heavier objects constituting the CUO could be white dwarfs, neutron stars, stellar black holes, or unresolved binaries.

We compared the radial distribution from our cleaned merged sample and from X-ray binaries from Bahramian et al. (2020). Figure 17 shows the cumulative distribution function (CDF) of these two datasets in the range of 2.7 − 100″, where our HST data was complete. It is clear that these two populations do not follow the same radial distribution (we find a KS p-value of 1.7 × 10⁻⁴). Furthermore, the bulk of the X-ray binaries seems to be located within 6″ − 50″ arcsec (Fig. 17). This is consistent with the CUO effective radius of 2.5 to 5″ (Table 3).

One may ask which among white dwarfs, neutron stars, BHs and massive binaries dominates the mass of the CUO. We can first discard unresolved binaries of main-sequence stars. Indeed, if they are unobserved, their total magnitude must be fainter than F606W = 22, which is the rough magnitude limit of our HST sample (Gaia is not relevant given the very low effective radius of the CUO). There is no way that their total mass can exceed that of the bulk of our sample. We could then have unresolved binaries of a main-sequence star with a compact star (white dwarf or neutron star) or possibly a BH. But the main-sequence star will have a mass of m_faint = 0.25 M_⊙ (Sect. 7.2.2), thus at least three times lower than our cut at m_bright = 0.77 M_⊙ and at least 20 times lower than that of a BH. Therefore, the main-sequence mass can be neglected and we are left with the compact star or BH.

We can compare the total mass in white dwarfs, neutron stars and BHs, by integrating over the zero-age (stellar mass function of the) main sequence (ZAMS):

$$\begin{array}{c}\hfill {\displaystyle {\int }_{m_{\min }}^{m_{\max }}{m}_{\mathrm{remnant}}(m)n(m)\mathrm{d}m,}\end{array}$$(29)

where n(m) is the ZAMS. We adopted the initial – final (remnant) mass relations from Eqs. (4)–(6) of Cummings et al. (2018) for white dwarfs and from Eqs. (C1), (C2), (C11), and (C15) of Spera et al. (2015) for neutron stars and BHs. Figure 18 displays these initial - final mass relations. We took a minimum remnant mass of 0.77 M_⊙, corresponding to the maximum main-sequence mass (Sect. 7.2.2), below which mass segregation is not effective against the most massive main-sequence stars. The maximum possible neutron star mass is M ≃ 2.15 M_⊙ (Rezzolla et al. 2018). We conservatively assumed a mass gap between 2.15 and 5 M_⊙ from the lack of LIGO detections in this mass range. We also considered a maximum stellar BH mass above which pair-instability supernovae fully explode the progenitor star leaving no remnant, adopting M_BH, max = 45 M_⊙ (Farmer et al. 2019) or 52 M_⊙ (Woosley 2017). We ignored the formation of very massive BHs above the pair-instability gap (i.e., M_BH >  133 M_⊙ Woosley 2017).

Summing the masses of each component by integrating Eq. (29) over the ZAMS mass function (i.e., the initial mass function, which in this mass range always has the Salpeter 1955 slope of −2.3) leads to BHs accounting for ≈58% of the CUO, with only ≈ 30% from white dwarfs and ≈ 12% from neutron stars. These fractions vary little with the maximum allowed BH mass: with 55% of the CUO mass in BHs with M_BH <  45 M_⊙ (Farmer et al. 2019) or 60% with M_BH <  52 M_⊙ (Woosley 2017). These fractions are insensitive to the metallicity of NGC 6397 from Z = 0.00013 to 0.0004 (i.e., from 1% to 3% solar, consistent with the metallicities derived for NGC 6397 by Marín-Franch et al. 2009 and Jain et al. 2020, respectively).

However, we must take this CUO mass dominance by BHs with caution because BHs may merge, losing of order of 5% of their mass to gravitational waves (e.g., Abbott et al. 2016), leading to kicks, some of which are strong enough to drive them out of the GC (Sect. 10.5). Still, one part of the BH population may merge and end up escaping the GC, while another part has not merged, but would nevertheless be located at low radii thanks to orbital decay by dynamical friction. If the CUO mass fraction in BHs were at least f₀ = 0.55 or 0.6 before merging and escapes, and if a mass fraction f_d of this BH component disappeared through mergers or escapes, then BHs would still dominate the CUO if f₀ (1 − f_d) > 1 − f₀, i.e. f_d <  2 − 1/f₀ = 0.19 or 0.33 considering the maximum BH mass of 45 M_⊙ (Farmer et al. 2019) or 52 M_⊙ (Woosley 2017). Put another way, BHs could contribute to half the CUO mass, if the maximum surviving BH mass is 40 M_⊙ according to our integrations of Eq. (29). The more massive BHs would sink faster to the center by dynamical friction and preferentially merge, leaving these lower mass ones. This discussion is a simplification because orbital decay by dynamical friction is stochastic, and one needs to test this with N-body simulations.

Our preference for the CUO model agrees with the suggestions by Zocchi et al. (2019) and Mann et al. (2019) that such a CUO can mimic an IMBH. Both studies obtained constraints on the mass of the stellar BH population in ω Cen and 47 Tuc, respectively. While those GCs are not core-collapsed, contrary to NGC 6397, it is interesting to compare our results with theirs. We determine that the CUO accounts for 1 to 2% of the GC mass, and if BHs do not merge or escape, the bulk of the CUO should be in BHs. Zocchi et al. (2019) tried different BH mass fractions, and their best fits were obtained with 1 to 5% of the GC mass in BHs. Mann et al. (2019) used as input a BH component whose mass is 1.4% of the GC mass, but they show that only 8.5% of this mass can be retained to avoid negative IMBH mass. Zocchi et al. (2019) did not provide a value for the scale radius of their stellar BHs, while Mann et al. (2019) provided a radius as an input parameter, which amounts to 6% of their fitted GC scale radius. In comparison, we are able to strongly constrain the CUO scale effective to 1.7 ± 0.5% of the GC effective radius, making our BH system much more concentrated.

A more robust conclusion of our integration of Eq. (29) is that white dwarfs always dominate the neutron stars, by a factor of ≈4, regardless of the possible dominance of BHs in the mass of the CUO component. Moreover, neutron stars can also merge together (or with black holes), and such an event has been detected by LIGO (Abbott et al. 2017). Presumably, the gravitational waves emitted will be much weaker and the lower momentum vector of the waves should lead to smaller kicks on the remnant (despite the lower masses of neutron stars compared to BHs). In summary, the CUO mass should be dominated by BHs, but must also contain white dwarfs, contributing four times more mass to the CUO than neutron stars.

Interestingly, stellar black holes in a dense inner subcluster, such as our CUO, could represent a major channel for the gravitational wave detections by LIGO/VIRGO (Portegies Zwart & McMillan 2000). These detections involve black holes of mass ≳20 M_⊙ (see Abbott et al. 2019) where electromagnetic wave detections are still lacking (e.g., Casares 2007). Black hole mergers should be even more frequent, presently, among the core-collapsed GCs such as NGC 6397, whose high inner densities allow for faster dynamical processes, such as dynamical friction, dynamical creation of black hole binaries (Heggie & Hut 2003), and subsequent hardening in three-body encounters (Heggie 1975). Our dynamical analysis provides a promising route to determine the locations of these stellar mass black holes.

10.7. Final thoughts

The future of GC and IMBH science is very exciting, thanks to Gaia now supplementing HST PMs at large projected radii. Our discovery of a diffuse central mass, composed in large part of stellar mass black holes, enrichens the physics of the inner regions of GCs, and renders the search of IMBHs in Milky Way GCs even more delicate. Continued pointings of GCs with HST and soon James Webb Space Telescope will lead to longer baselines and more accurate PMs. The third data release of the Gaia mission will double the PM precision, thus enabling more accurate mass-orbit modeling, not only of nearby GCs such as NGC 6397, but also of more distant ones, in conjunction with HST data. Future gravitational wave missions such as the Laser Interferometer Space Antenna (LISA) will probe IMBHs as well as stellar BHs, including in our Milky Way (Sesana et al. 2020). The physics of the inner parts of GCs with possible IMBHs as well as subclusters of stellar-mass black holes is truly enticing.

---

Acknowledgments

We are most grateful to Andrea Bellini and Sebastian Kamann for respectively providing PM and LOS velocity data for NGC 6397, without which this paper would not have been possible. We also thank Thomas de Boer and Mark Gieles for providing their unpublished estimate of the projected half-light radius of NGC 6397. We acknowledge the anonymous referee who provided useful comments and references. We also thank Gwenaël Boué and Radek Wojtak for useful comments and suggestions in our mass modeling approach, Lukas Furtak for his help with the PARSEC code, Marta Volonteri for insightful discussions, and Pierre Boldrini, Laura Watkins and Natalie Webb for useful references. Eduardo Vitral was funded by a grant from the Centre National d’Études Spatiales (CNES) and an AMX doctoral grant from École Polytechnique. This work greatly benefited from public software packages SCIPY, ASTROPY (Astropy Collaboration 2013), NUMPY (van der Walt et al. 2011), COSMOMC (Lewis & Bridle 2002), EMCEE (Goodman & Weare 2010), MATPLOTLIB (Hunter 2007), as well as Pierre Raybaut for developing the SPYDER Integrated Development Environment. We also intensively used TOPCAT (Taylor 2005).

References

Appendix A: Deprojection of the Sérsic surface density profile

The Sérsic surface density model used in this work is deprojected in the same way as Vitral & Mamon (2020, hereafter VM20), but we extend the deprojection to lower radii: 10⁻⁴ R_e instead of 10⁻³ R_e. We consider again the approximations of Lima Neto et al. (1999, hereafter LGM99), Simonneau & Prada (2004, hereafter SP04), as well as the formula of VM20, to find which is better suited for each region of the [n × r/R_e] domain.

A.1. New coefficients for VM20 approximation extending to very low radii

The coefficients of the VM20 approximation were originally calculated for a logarithmic grid of [n × r/R_e] with −3 ≤ log(r/R_e)≤3 (100 steps) and 0.5 ≤ n ≤ 10 (50 steps). However, our HST data could extend to even lower radii than 0.001 R_e, if R_e >  7′ (the lowest projected radius is $0\stackrel{″}{.}42$). Indeed, such large effective radii were found for the faint stars of the two-population models 14 and 15. We therefore recomputed the VM20 approximation, using a different region of the [n × r/R_e] domain, and an even finer grid.

We used the same approach as we did in VM20, but this time we performed the numerical deprojection of the Sérsic profile with MATHEMATICA 12, instead of PYTHON. Given Fig. 4 of VM20 and the lower limit of r/R_e we needed to attain, we calculated the best VM20 parameters for a new logarithmic spaced region [n × r/R_e] limited to −4 ≤ log(r/R_e)≤3 (150 steps) and 0.5 ≤ n ≤ 3.5 (100 steps). The resulting coefficients a_i, j, presented in the same way as VM20 (see their Eq. (28)), can be found online¹⁹ and we hereafter refer to them as VM20bis.

A.2. Choice of deprojection approximation

The two left panels of Fig. A.1 display the best fitting approximations in 100 × 150 grid in [log n × log (r/R_e)]. We included in MAMPOSST-PM a simplified choice of best approximations to the deprojected Sérsic mass and density profiles. For the mass profile we used

• n × r/R_e ∈ [0.5; 1.5] × [1; 10³]: LGM99

• n × r/R_e ∈ [0.5; 3.4] × [10⁻⁴; 1): VM20bis

• n × r/R_e ∈ (3.4; 10] × [10⁻⁴; 1)∪(1.5; 10]×[1; 10³]: SP04

For the density profile, the division was even simpler:

• n × r/R_e ∈ (3.4; 10] × [10⁻⁴; 10³]: SP04

• n × r/R_e ∈ [0.5; 3.4] × [10⁻⁴; 10³]: VM20bis

The outcome of this approximation is highly accurate, as seen in the right panels of Fig. A.1 (which are the equivalent of Fig. 3 in VM20), with the grid lines indicating the divisions adopted for the 3 approximations. The reader can verify that the relative precision of this hybrid model is on the order of ∼0.1%.

Appendix B: Handling of field stars

B.1. Proper motion distribution function

In this appendix, we justify the need for using wider tails than a Gaussian for the distribution of PMs in the field stars, as mentioned in Sect. 5.2.3.

We downloaded Gaia DR2 PM data in four regions around NGC 6397 (and also for two other GCs, M 4 and NGC 6752, in order to check for generality). These four regions were chosen by doing a cone search, with a 30′ aperture, for positions 5° distant from the GC centers (α_GC, δ_GC), north, south, east, and west. We applied the same quality flags that we had applied to the NGC 6397 data (Sect. 5.2.1): the inequalities of Eqs. (12), (16) and (17).

We then fitted the distribution of PM moduli using both a Gaussian and the form of Eq. (18). We estimated the mean μ_α, * and mean μ_δ for both distribution functions, a dispersion σ for the Gaussian assumption and for Eq. (18), we estimated the scale radius and outer slope. Finally, we took into account the convolution of both distributions with Gaussian errors.

Figure B.1 shows the distribution of PM moduli, for the four regions NGC 6397. One can easily verify that the new expression fits extremely better than a Gaussian in Fig. B.1. Moreover, the calculated kurtosis excess of both μ_α, * and μ_δ always gave huge values (from 16 to over 400, compared to 0 expected for a Gaussian), which clearly implies a non-Gaussian behavior.

To check for robustness of our fits, we also verified that whenever limiting the fitting range to exclude the wider wings of the FS distribution (which cannot be done when considering a GC relatively separated from the FS in PM space), a Gaussian distribution is well fitted. Thus, we decided to employ Eq. (18) throughout our study in order to account for the entire PM range.

B.2. Convolution of field stars distribution

For the GC component, the local Gaussian-like velocity distribution function (VDF) must be convolved by the Gaussian LOS and POS velocity errors, as in Eq. (7). Similarly, since the field star LOS VDF is a Gaussian, the LOS velocity errors are added in quadrature to the LOS velocity dispersion.

On the other hand, since the interloper PM distribution function (Eq. (18)) is not Gaussian-like shaped, we need to perform the integral of its convolution. Calling R and R_o the respective true and observed PMs and ϵ the PM error, it is straightforward to follow the recipe of Binney & Mamon (1982, Appendix C) for the convolution of two-dimensional data with circular symmetry:

$$\begin{array}{c}\hfill p_{\mathrm{conv}}(R_{\mathrm{o}})=\frac{R_{\mathrm{o}}}{{ϵ}^2}{\displaystyle {\int }_0^{\infty }\mathrm{d}Rp(R)exp(-\frac{R^2+R_{\mathrm{o}}^2}{2{ϵ}^2})I_0\left(\frac{R{R}_{\mathrm{o}}}{{ϵ}^2}\right)·}\end{array}$$(B.1)

Mamon & Vitral (in prep.) provide an approximation for this integral.

All Tables

All Figures