Hubble Space Telescope Proper Motion (HSTPROMO) Catalogs of Galactic Globular Clusters. V. The Rapid Rotation of 47 Tuc Traced and Modeled in Three Dimensions^*

Globular cluster (GC) formation, internal dynamical evolution, and the effects of the external tidal field of the host galaxy are expected to leave a number of fingerprints on a cluster's structural, morphological, and kinematical properties.

Until recently, the dynamical characterization of most GCs has been limited to the surface brightness or projected star count radial profiles. These profiles can shed light on some aspects of cluster dynamical evolution—such as the identification of the systems that have already evolved past the core-collapse phase (see, e.g., Djorgovski & King 1986; Chernoff & Djorgovski 1989; Djorgovski & Meylan 1994; Trager et al. 1995), and the calculation of a number of fundamental global properties (see, e.g., McLaughlin & van der Marel 2005; Miocchi et al. 2013)—but they provide only a very partial view of a cluster's dynamical state. To construct a complete dynamical picture of clusters, we require an accurate understanding of both their internal kinematics and their detailed morphologies. Such information places key constraints on theoretical studies and thus allows us to reconstruct cluster dynamical histories and determine the role played by different dynamical processes (e.g., Giersz & Heggie 2011).

After some early pioneering work on cluster kinematical properties (see, e.g., Pryor & Meylan 1993; Meylan & Heggie 1997, and references therein; van Leeuwen et al. 2000), there has been a recent revival in observational studies of internal cluster kinematics based on ESO/VLT radial velocities (see, e.g., Bellazzini et al. 2012; Lanzoni et al. 2013; Kamann et al. 2016) and Hubble Space Telescope (HST) based proper-motion (PM) measurements (see, e.g., Bellini et al. 2014; Watkins et al. 2015a, hereafter Paper I and Paper II, respectively).

These kinematics have already been used to address a number of fundamental issues, including the existence of intermediate-mass black holes at the centers of GCs (see, e.g., Anderson & van der Marel 2010; Noyola et al. 2010; van der Marel & Anderson 2010; Lanzoni et al. 2013; Lützgendorf et al. 2013), the kinematical differences between multiple stellar populations (Richer et al. 2013; Bellini et al. 2015; Cordero et al. 2017), differences between the velocity dispersion of stars with different masses and the degree of energy equipartition (see, e.g., Baldwin et al. 2016, Paper IV), the presence and origin of anisotropy in the velocity distribution (Paper II), and the strength of cluster rotation and its possible link to the cluster morphology (see, e.g., Bellazzini et al. 2012; Bianchini et al. 2013; Kacharov et al. 2014; Fabricius et al. 2014; Lardo et al. 2015; Boberg et al. 2017). Forthcoming data from the Gaia mission (see, e.g., Pancino et al. 2017) will further enrich the observational landscape.

On the theoretical side, renewed efforts are being made to expand the numerical and analytical tools necessary to interpret the results of these observational studies and gather a deeper understanding of their implications in the context of GC dynamical evolution. These efforts include studies proposing new distribution-function-based models of rotating and anisotropic models (see, e.g., Varri & Bertin 2012; Gieles & Zocchi 2015; de Vita et al. 2016), extending the widely used (but limited to spherical symmetry and isotropic velocity distribution) models such as the King (1966) or Wilson (1975) models, and numerous numerical studies exploring, for example, the evolution of rotating and anisotropic models (see, e.g., the early studies of Einsel & Spurzem 1999; Kim et al. 2002; Ernst et al. 2007; and the more recent investigations by Hong et al. 2013; Tiongco et al. 2016b, 2017; Zocchi et al. 2016), the evolution of the degree of energy equipartition (see, e.g., Trenti & van der Marel 2013; Bianchini et al. 2016; Webb & Vesperini 2017), the role of potential escapers in the observed velocity dispersion profiles (see, e.g., Küpper et al. 2010; Claydon et al. 2017), and the kinematical implications of star loss from GCs (see, e.g., Keenan & Innanen 1975; Tiongco et al. 2016a).

On the observational, PM-based side, measurements of GC rotation in the plane of the sky had been completely lacking until van Leeuwen et al. (2000) took advantage of a large number of photographic plates spanning half a century to measure a plane-of-the-sky differential rotation for the GC ω Centauri (see their Figure 18). A few years later, Anderson & King (2003a) took advantage of the high-precision astrometric capabilities of HST to obtain a direct measurement of the plane-of-the-sky rotation of 47 Tuc (NGC 104) using two diametrically opposite fields ∼5' from the cluster center. The authors exploited the unique feature of 47 Tuc to have a large number of relatively bright Small Magellanic Cloud (SMC) stars in the background, and SMC stars themselves could be used as reference objects. In contrast, much fainter galaxies or much rarer QSOs must be used as reference objects for the vast majority of the other GCs. A few other GCs analogous to 47 Tuc are (1) NGC 362, which is also in front of SMC stars and will be the subject of a forthcoming paper in this series; (2) NGC 6652, which is in front of Sagittarius dwarf spheroidal (Sgr dSph) galaxy stars; and (3) NGC 6681, also in front of Sgr dSph stars. The plane-of-the-sky rotation of NGC 6681 has indeed been recently measured (with a null result) by Massari et al. (2013), but using background galaxies as a reference. The Massari et al. (2013) work represents the third—and so far the last—observational, PM-based work on the subject.

In this paper, we extend the pioneering work of Anderson & King (2003a) and measure, for the first time, the rotation curve of 47 Tuc in the plane of the sky from the center of the cluster out to about 13' (about four half-light radii; Harris 1996). We present state-of-the-art, high-precision, HST-based PM measurements for the cluster, which we combine with existing line-of-sight (LOS) velocities to determine the strength of the cluster present-day rotation and anisotropy in the velocity distribution. By combining the kinematical data with the known surface brightness profile, we build a detailed self-consistent model of the cluster. This work represents the most complete dynamical characterization of any Galactic GC and clearly illustrates how only such a detailed study can reveal the intrinsic dynamical properties of a cluster and provides the needed constraints to explore the possible evolutionary paths leading to the present-day observed properties.

The structure of this paper is as follows: In Section 2 we present the HST data set and the reduction techniques used to measure high-precision PMs. The results concerning the PM-based rotation of 47 Tuc in the plane of the sky and the velocity anisotropy are presented, respectively, in Sections 3 and 4. In Section 5 we present the results of a detailed dynamical model fitting to the observational data using the distribution-function-based models of Varri & Bertin (2012). Conclusions are summarized in Section 6.

We measured PMs in four fields located at different radial distances from the center of 47 Tuc. These fields are roughly aligned along the same axis with respect to the cluster's center and span over 1200'' (>20') on the sky. In Figure 1 we show the footprints of these four fields superimposed on a Sloan Digital Sky Survey (SDSS) image of the cluster. The four fields are hereafter identified, from east to west, as the inner field, the central field, the calibration field, and the outer field. We made use of exposures taken from three different HST cameras: the Wide-Field Planetary Camera 2 (WFPC2), both the High-Resolution Channel (HRC) and the Wide-Field Channel (WFC) of the Advanced Camera for Surveys (ACS), and the Ultraviolet-VISible (UVIS) channel of the Wide-Field Camera 3 (WFC3).⁹

Zoom In Zoom Out Reset image size

The data reduction and analysis are based on _flt-type exposures (or the equivalent _c0f format for WFPC2), as they preserve the un-resampled pixel data for optimal stellar profile fitting. All ACS/WFC and WFC3/UVIS exposures were corrected for charge transfer efficiency (CTE) defects (Anderson & Bedin 2010).

Images were then reduced using the software family img2xym (Anderson & King 2006a), employing either single (HRC), spatially varying (WFPC2), or spatially and time-dependent empirical point-spread functions (WFC and UVIS; Anderson & King 2000, 2006a, 2006b; Bellini et al. 2013). Stellar positions were corrected for geometric distortion using the state-of-the-art solutions provided by Anderson & King (2003b, 2006a, 2006b), Bellini & Bedin (2009), and Bellini et al. (2011). The photometric data of the central field come directly from Anderson et al. (2008). Photometry of the other fields was calibrated following the prescriptions given in Sirianni et al. (2005).

Stellar positions in each exposure were transformed into a common, distortion-free reference frame (the master frame). Only stars measured in at least four distinct exposures with a time baseline of at least 6 months were considered for the present analysis. PMs were computed using the central overlap method (Eichhorn & Jefferys 1971), in which each exposure is considered as a stand-alone epoch. In a nutshell, we transformed stellar positions as measured on the individual exposures on the master frame, by means of general, six-parameter linear transformations. For each star, its master-frame-transformed positions as a function of the epoch are fitted by a straight line, the slope of which gives us a direct measurement of the stellar motion. We applied a careful data rejection procedure to remove outliers or mismatches (see Section 5.5 of Paper I for more details). Extensive simulations have demonstrated the reliability of both estimated PMs and PM errors.

In the following, we describe the reduction procedures of each field.

2.1. The Central Field

2.2. The Calibration Field

This particular field was selected as one of the instrument calibration fields for the ACS and WFC3 detectors and has been repeatedly observed since 2002, when ACS was installed on board HST. We chose a subsample of deep exposures taken from 2002 to 2014 with ACS/WFC and WFC3/UVIS (see Table 1). The final PM catalog contains 13,466 objects, 12,132 of which passed our high-quality selection criteria.

Table 1.  List of Observations Used for the Calibration Field

Epoch	GO	Instr.	Filter	$N\times \exp .\mathrm{time}$
2002.30	9018	ACS/WFC	F606W	1 × 765 s, 1 × 1200 s
			F814W	1 × 690 s, 1 × 1020 s
2006.67	10730	ACS/WFC	F435W	8 × 350 s
	10737	ACS/WFC	F435W	1 × 339 s
			F555W	1 × 339 s
			F606W	1 × 339 s
			F814W	1 × 339 s
2009.54	11444	WFC3/UVIS	F606W	24 × 350 s
	11452	WFC3/UVIS	F814W	6 × 350 s
2012.22	12692	WFC3/UVIS	F606W	2 × 350 s
2014.35	13596	ACS/WFC	F435W	2 × 339 s
			F475W	2 × 339 s
			F606W	2 × 339 s
			F775W	2 × 339 s
			F814W	2 × 339 s
			F850L	2 × 339 s

Download table as:  ASCIITypeset image

2.3. Inner and Outer Fields

The inner and outer fields were observed with either ACS/WFC or WFC3/UVIS in only one epoch. All other available observations were taken with the WFPC2. The inner field consists of one single ACS pointing and all the available WFPC2 exposures that overlap it. The outer field consists of three marginally overlapping UVIS pointings, roughly aligned along the decl. direction, and all the available WFPC2 exposures overlapping them. The complete lists of observations for the inner and outer fields are reported in Tables 2 and 3, respectively.

Table 2.  List of Observations Used for the Inner Field

Epoch	GO	Instr.	Filter	$N\times \exp .\mathrm{time}$
1994.75	5370	WFPC2	F606W	4 × 600 s, 7 × 1000 s
			F814W	7 × 1000 s
1994.86	5369	WFPC2	F606W	9 × 300 s, 1 × 500 s
1999.79	8095	WFPC2	F606W	1 × 400 s, 7 × 500 s
2013.34	12971	ACS/WFC	F435W	1 × 290 s, 1 × 690 s
			F555W	1 × 360 s, 1 × 660 s

Download table as:  ASCIITypeset image

Table 3.  List of Observations of the Outer Field

Epoch	GO	Instr.	Filter	N × exp.time
2002.35	9318	WFPC2	F606W	1 × 300 s, 1 × 400 s
				7 × 500 s
2002.40	8059	WFPC2	F450W	1 × 300 s
			F606W	1 × 160 s
			F814W	1 × 100 s, 1 × 300 s
2003.40	9709	WFPC2	F606W	1 × 400 s, 7 × 500 s
2010.34	11677	WFC3/UVIS	F390W	1 × 1048 s, 1 × 1099 s
				1 × 1212 s, 1 × 1355 s
				2 × 1400 s
			F606W	1 × 1252 s, 3 × 1347 s
				1 × 1398 s, 1 × 1402 s

Download table as:  ASCIITypeset image

The PM reduction tools presented in Paper I were not designed to be used on WFPC2 data, because the astrometric precision reachable with WFPC2 exposures is generally far less than what we can achieve with later HST optical imagers. The WFPC2 camera was composed of four 800 × 800 pixel 12-bit detectors, with a pixel scale of 99.6 mas pixel⁻¹ for the three WF chips and 45.5 mas pixel⁻¹ for the PC chip. All four chips suffered from ∼16% vignetting. The incident light was split into four beams (one per chip) by a four-faced pyramid mirror, which provided a challenging registration of the observed relative chip positions. Finally, there is no CTE correction available for the WFPC2.

Nevertheless, WFPC2 observations represent the only available first epochs for the inner and outer fields and allowed us to compute PMs with a time baseline of 18.6 and 8 yr, respectively. Our PM reduction tools are scalable, and it is straightforward to include data coming from different instruments/cameras. We treated each WFPC2 chip independently, so as to minimize interchip transformation errors. We modeled single-exposure expected astrometric errors as a function of the instrumental magnitude using all the available WFPC2 exposures in the two fields. The expected errors are used as a first-guess weight during the PM-fitting procedures (see Section 5.2 of Paper I for more details).

The computed PMs for the inner and the outer fields are in agreement with those computed in the central and calibration fields in terms of expected intrinsic values and associated errors. On the other hand, the available data set does not allow us to adequately study and minimize the impact of systematic effects in the quoted PM errors, as we did for the ACS and UVIS detectors. As a result, we cannot as reliably study the internal kinematics of 47 Tuc stars (which relies on the subtraction in quadrature of PM errors) using inner- and outer-field measurements as for the central and calibration fields.

The inner-field PM catalog contains 2187 objects, 2084 of which passed our high-quality selection criteria. The outer field consists of 648 total objects, of which 579 are identified as high-quality measurements.

The top panels of Figure 2 show the field of view of each of the four fields around the center of the cluster. From left to right we have (1) the inner field, (2) the central field, (3) the calibration field, and (4) the outer field. Concentric circles, in red, give an idea of the radial extension of the data.

Zoom In Zoom Out Reset image size

HST exposures all have different roll angles for different epochs, and axis rotation is one of the six parameters that are solved for when we transform stellar positions as measured on single exposures into the master frame. Because of this, any direct sign of cluster rotation (if present) is absorbed by the linear terms of the coordinate transformations. Since our PMs are relative to the bulk motion of the cluster, any other object in the field that is not a cluster member would have a systematic component in its PM measurement that is equal in size and with the opposite sign of the bulk rotation of the cluster.

This systematic PM component is clearly visible in the PM diagrams shown in the middle panel of Figure 2. Cluster stars are in black, while background SMC stars are marked as red crosses. As a reference, a black crosshair in each panel highlights the barycenter of SMC stars as measured in the central field, where, because of symmetry, the systematic PM rotation component is minimized. The white crosses near the center of the crosshairs mark the median loci of SMC stars in each of the four fields. The mean PM of SMC stars in the inner field is shifted toward smaller ${\mu }_{\delta }$ values, while that of the calibration and outer fields is shifted toward larger values of ${\mu }_{\delta }$. It is clear from the figure that SMC stars appear to rotate counterclockwise with respect to the center of 47 Tuc, which directly translates into a clockwise rotation of 47 Tuc in the plane of the sky.

The bottom panels of Figure 2 show the color–magnitude diagrams (CMDs) of 47 Tuc stars (black dots) and SMC stars (red crosses) in different magnitude/color combinations for the four fields. To select bona fide SMC stars, we took advantage of their location on both the PM diagrams and the CMDs. SMC stars on the PM diagram are clustered around a well-defined location, clearly distinct from that occupied by 47 Tuc stars. We preliminary selected SMC candidates as all the sources 4σ outside the distribution of 47 Tuc stars. This rough selection necessarily includes field stars and a few cluster members. We then iteratively computed 4σ-clipped median values ${\mu }_{\alpha }\cos \delta$ and ${\mu }_{\delta }$ of the barycenter of SMC stars in the PM diagram and kept only stars within 4σ from this location. Finally, we excluded by eye from the SMC selection all those stars that did not lie in the SMC region on the CMDs. We repeated one more time the process of computing the location of the barycenter of SMC stars in the PM diagram, and we kept only SMC stars within 4σ from their median location.

The final sample of SMC stars is what we show in red in Figure 2. We have 96 SMC stars in the inner field, 164 in the central field, 1495 in the calibration field, and 235 in the outer field.

In order to properly measure the apparent rotation of SMC stars, we need to define a reference point (the zero-point) on the PM diagram with respect to which to compute their tangential component of the motion (${\mu }_{\mathrm{TAN}}$). To do this, we selected all SMC stars within the largest circle that can be fully contained within the central field and determined a 3σ-clipped estimate of their barycenter location on the PM diagram. This location has coordinates ${\mu }_{\alpha }\cos \delta =-4.797\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$, ${\mu }_{\delta }=1.244\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$, and it defines the center of the crosshairs in the middle panels of Figure 2.

We zero-pointed the PM of SMC stars to the reference point defined above and computed their average ${\mu }_{\mathrm{TAN}}$ at different equally populated radial intervals with respect to the cluster's center. A finer radial subdivision (about 28 stars per bin) is applied to the central field, where the fastest rotational variation is found (at the cost of larger measurement errors). We then derived a single rotation measurement for the inner field (96 stars), 10 measurements for the calibration field (about 150 stars each), and two measurements for the outer field (117 stars each). The computed quantities, as a function of the radial distance, are shown with error bars in Figure 3 and are listed in Table 4. We adopted the convention of using negative ${\mu }_{\mathrm{TAN}}$ values in case of a clockwise rotation in the plane of the sky (north to west). The horizontal bars in Figure 3 indicate the size of the radial intervals within which each ${\mu }_{\mathrm{TAN}}$ value is computed. For completeness, we included (in red) the rotation value computed by Anderson & King (2003a). All these quantities are listed in Table 4. The cluster has a solid-body-like rotation within our central field. A simple least-squares fit of the form ${\mu }_{\mathrm{TAN}}={mr}$ gives m = −1.65 ± 0.10 $\mathrm{mas}\,{\mathrm{yr}}^{-1}$ arcsec⁻¹. Then, the rotation profile slowly flattens at larger radii, to reach the highest value of −0.312 ± 0.015 $\mathrm{mas}\,{\mathrm{yr}}^{-1}$ at ∼390'' (∼652) from the cluster center. The measured rotation slows down to about −0.21 $\mathrm{mas}\,{\mathrm{yr}}^{-1}$ in the outermost regions probed by our data.

Zoom In Zoom Out Reset image size

Table 4.  Rotation of 47 Tuc in the Plane of the Sky

$\langle r\rangle$	N${}_{\mathrm{STAR}}$	$\langle {\mu }_{\mathrm{TAN}}\rangle$	$\langle {\sigma }_{{\mu }_{\mathrm{TAN}}}\rangle$	$\langle r\rangle$	N${}_{\mathrm{STAR}}$	$\langle {\mu }_{\mathrm{TAN}}\rangle$	$\langle {\sigma }_{{\mu }_{\mathrm{TAN}}}\rangle$
('')		(mas yr)	($\mathrm{mas}\,{\mathrm{yr}}^{-1}$)	(arcsec)		($\mathrm{mas}\,{\mathrm{yr}}^{-1}$)	($\mathrm{mas}\,{\mathrm{yr}}^{-1}$)
62.39	28	−0.112	0.050	373.12	150	−0.296	0.013
83.97	28	−0.132	0.046	391.51	150	−0.312	0.015
94.92	28	−0.132	0.042	408.38	150	−0.301	0.013
101.42	28	−0.142	0.041	424.24	150	−0.288	0.017
109.86	28	−0.210	0.045	440.00	150	−0.260	0.013
124.28	24	−0.231	0.051	455.87	150	−0.271	0.014
219.02	96	−0.252	0.031	472.92	150	−0.246	0.014
321.49	150	−0.288	0.015	494.26	145	−0.275	0.013
(342)a		(−0.233)a	(0.055)a	705.53	118	−0.232	0.041
352.61	150	−0.283	0.016	775.05	117	−0.208	0.045

Note.

^aFrom Anderson & King (2003a).

Download table as:  ASCIITypeset image

The general shape of the rotation curve of 47 Tuc we have measured on the plane of the sky is qualitatively similar to that obtained by Bianchini et al. (2013) using LOS measurements. There is, though, an important difference between the two rotation profiles: the PM-based profile has a rotation peak in the plane of the sky that is about twice as large as that measured with LOS velocities, which we model in Section 5.

To give the reader a better sense of how the cluster is rotating on the plane of the sky, we show in Figure 4 the rotation map of the cluster derived with the data sets we analyzed. To obtain the map, we divided the data sets into either four quadrants (for the central field) or equally populated regions (for the other three fields) and computed the zero-pointed median ${\mu }_{\alpha }\cos \delta$, ${\mu }_{\delta }$ components of the motion of SMC stars in each region. These values (with the opposite sign) are shown as vectors departing from the median location of SMC stars within each region. The length of the vectors, in $\mathrm{mas}\,{\mathrm{yr}}^{-1}$, has been magnified by a factor of 750, for clarity. For completeness, we show in red the rotation axis as measured from LOS data (136° north to east; Bianchini et al. 2013).

Zoom In Zoom Out Reset image size

The careful reader might have noticed that the distribution of 47 Tuc stars in the PM diagram of the central field is rather circular, while this distribution is more flattened in the inner and calibration fields, to become again somewhat more circular in the outer field (middle panels of Figure 2). This is the effect of velocity dispersion anisotropy.

To properly measure the degree of velocity dispersion anisotropy of the cluster as a function of radius, we proceeded as follows. First, we need to select stars with reliable PMs and similar masses in the four fields. As we mentioned earlier, PMs computed in the inner and outer fields are based on WFPC2 measurements and might be affected by significant systematic effects. Our data sets do not allow us to adequately study and minimize the impact of systematic effects due to WFPC2 measurements. While these systematic errors are expected to only have second-order effects in the quoted PMs, they can significantly alter the estimated PM errors. Since PM errors are subtracted in quadrature when we want to compute velocity dispersion profiles, under/overestimating PM errors could lead to incorrect profiles. In what follows, we will include measurements coming from the inner and the outer fields only for completeness.

To select stars of similar mass in the four fields, we limit our selections to the magnitude range $20\lt {m}_{{\rm{F}}606{\rm{W}}}\lt 22$. The radial and tangential velocity dispersion profiles were estimated using the same method as in van der Marel & Anderson (2010), which corrects the observed scatter for the individual PM uncertainties.

The top panel of Figure 5 shows the radial velocity dispersion profile ${\sigma }_{\mathrm{RAD}}$ of 47 Tuc as a function of radius. In red we report the values computed in the inner and outer fields, for completeness. The horizontal bars illustrate the radial extent over which each point is obtained. The central radial velocity dispersion is about 0.64 $\mathrm{mas}\,{\mathrm{yr}}^{-1}$, or 13.65 $\mathrm{km}\,{{\rm{s}}}^{-1}$ at a distance of 4.5 kpc (to be compared to 11.5 $\mathrm{km}\,{{\rm{s}}}^{-1}$ for evolved stars; Harris 1996).¹⁰ The profile reaches a minimum value of about 0.31 $\mathrm{mas}\,{\mathrm{yr}}^{-1}$ in the outermost point of the calibration field (∼825, or about 495''). The middle panel of the figure shows the tangential velocity dispersion profile ${\sigma }_{\mathrm{TAN}}$ computed in the same radial intervals as for ${\sigma }_{\mathrm{RAD}}$. The two vertical lines mark the location of the core radius ${r}_{{\rm{c}}}=0\buildrel{\prime\prime}\over{.} 36$ (216) and the half-light radius ${r}_{{\rm{h}}}=3\buildrel{\,\prime}\over{.} 17$ (1902) (Harris 1996). Table 5 lists the values of the data points shown in Figure 5.

Zoom In Zoom Out Reset image size

Table 5.  Radial and Tangential Velocity Dispersion Profile Data of 47 Tuc

$\langle r\rangle$	${N}_{\mathrm{STAR}}$	$\langle {\sigma }_{\mathrm{RAD}}\rangle$	err${}_{{\sigma }_{\mathrm{RAD}}}$	$\langle {\sigma }_{\mathrm{TAN}}\rangle$	err${}_{{\sigma }_{\mathrm{TAN}}}$
(arcsec)		($\mathrm{mas}\,{\mathrm{yr}}^{-1}$)	($\mathrm{mas}\,{\mathrm{yr}}^{-1}$)	($\mathrm{mas}\,{\mathrm{yr}}^{-1}$)	($\mathrm{mas}\,{\mathrm{yr}}^{-1}$)
Central Field
31.15	210	13.39	0.64	13.68	0.66
37.15	209	14.16	0.74	14.94	0.78
41.45	308	13.22	0.54	14.58	0.59
46.53	461	13.84	0.45	13.95	0.45
51.64	540	13.38	0.44	13.15	0.43
56.52	675	13.55	0.35	13.32	0.34
66.95	2589	13.48	0.20	13.26	0.20
81.60	3305	13.11	0.17	12.54	0.15
96.49	3588	12.77	0.13	12.16	0.13
110.30	2160	12.19	0.18	11.93	0.18
123.46	671	11.81	0.35	11.74	0.35
134.10	143	12.40	0.87	11.53	0.81
Inner Field
178.30	220	10.04	0.49	9.11	0.44
214.73	220	9.59	0.50	8.88	0.47
263.52	213	9.56	0.42	7.72	0.34
Calibration Field
304.81	147	9.41	0.79	8.27	0.48
324.16	147	8.89	0.44	8.28	0.48
341.59	147	8.32	0.31	7.35	0.43
357.14	147	8.68	0.52	7.30	0.43
372.00	147	8.17	0.32	6.96	0.40
386.25	147	8.72	0.45	7.11	0.42
400.55	147	7.71	0.28	7.32	0.43
413.98	147	7.96	0.40	6.09	0.36
430.15	147	8.55	0.41	7.24	0.42
448.01	147	7.49	0.53	6.94	0.41
468.00	147	7.74	0.60	7.46	0.44
494.83	138	6.68	0.25	6.43	0.39
Outer Field
699.00	78	7.89	0.73	7.20	0.68
774.47	77	7.86	0.73	7.25	0.68

Download table as:  ASCIITypeset image

The deviation from isotropy (${\sigma }_{\mathrm{TAN}}$/${\sigma }_{\mathrm{RAD}}-1$) as a function of radius is shown in the bottom panel of Figure 5. The horizontal line at 0 indicates an isotropic system. The center of the cluster is isotropic, with increasing radial anisotropy moving outward. It is worth noting that this trend agrees with what we saw in Paper II for the entire sample of 22 GCs, albeit over a smaller radial range. The degree of anisotropy is significant in the calibration field, i.e., at radial distances between 5' (300'') and 8' (480'').

Both the radial and the tangential velocity dispersion profiles of 47 Tuc seem to drop in the centermost radial bin. The first and second points of the dispersion profiles are nonetheless consistent with each other to within less than 1σ. Changing the size of the first few radial bins produced similar trends, all consistent with being flat in the centermost regions within the error bars. Also note that these profiles are based on relatively faint main-sequence (MS) stars, which suffer the most from the highly crowded conditions of the core of the cluster. The velocity dispersion profile of 47 Tuc we published in Paper II, based on much brighter, higher signal-to-noise ratio RGB stars, does not show any central drop.

The goal of this section is to provide a dynamical model to describe the comprehensive set of available observations, comprising the PM kinematics analyzed above in addition to the classic LOS data and photometry. We will use a family of physically motivated distribution-function-based models (Varri & Bertin 2012), recently applied to a selected sample of Galactic GCs (Bianchini et al. 2013; Kacharov et al. 2014).

These self-consistent models have been specifically constructed to describe quasi-relaxed stellar systems with realistic differential rotation, axisymmetry, and pressure anisotropy. The models are defined by four dimensionless parameters (concentration parameter Ψ, rotation strength parameter χ, and the parameters b and c determining the shape of the rotation profile). A full description of the distribution function and of the parameter space is given in Varri & Bertin (2012). Since these models allow only for a single-mass component, they do not take into consideration the effects connected with mass segregation. The implications for this assumption are described in Section 5.1.1.

We will use the PM-based profiles described in Section 4 (tangential and radial velocity dispersion profiles and rotation profile in the plane of the sky), the LOS-based velocity profiles (velocity dispersion and rotation profiles) reported in Bianchini et al. (2013), and the photometric data from Trager et al. (1995). We do not consider in the fit procedure the data obtained from the inner and outer fields, as described in Section 4.

5.1. Fitting Procedure

We follow the same fitting procedure outlined in Kacharov et al. (2014) and Bianchini et al. (2013), as summarized in the following. The comparison between the differentially rotating axisymmetric models and the observations requires us to specify four dimensionless parameters and five additional quantities: three physical scales (i.e., the radial scale r₀, the central surface density ${{\rm{\Sigma }}}_{0}$, and the velocity scale v₀), the inclination angle i between the rotation axis and the LOS direction, and the cluster distance. Since the adopted dynamical models are characterized by deviations from isotropy in configuration and velocity space, the choice of the inclination angle plays a fundamental role in the fitting procedure. To exploit such an additional degree of freedom, we initially assume the value i = 30°, derived by applying Equation (8) of van de Ven et al. (2006) (linking the average motion along the LOS to the one in the plane of the sky) to our data. We will later explore how the best-fit parameters change with different inclinations (i = 25°–35°). Since PMs are measured in mas yr⁻¹ in the plane of the sky, a multiplying factor of 4.74d is needed to convert these values to $\mathrm{km}\,{{\rm{s}}}^{-1}$. We will assume a distance of 47 Tuc of $d=4.5\,\mathrm{kpc}$ (Harris 1996).

The fitting procedure is twofold. First, we determine the dimensionless parameters needed to reproduce the observed value of V/σ and the observed position of the rotation peak (for further details see Sections 3.1 and 3.5 of Bianchini et al. 2013). Then we calculate the physical scales through χ² minimization for all the kinematic profiles (three dispersion profiles, two rotation curves) to obtain the radial scale r₀ and the velocity scale v₀, and then we fit the surface density profile to obtain the central surface density Σ₀. This provides all the constraints needed to determine the best-fit dynamical model.

We wish to emphasize that, during the fitting procedure, the projection of the self-consistent dynamical models is performed by sampling from the relevant distribution function a discrete set of N = 2,048,000 particles and then performing a rotation of such a discrete system to match the relevant inclination angle. The theoretical kinematic and photometric profiles are then calculated by following the same procedures applied for the construction of the observational profiles (i.e., by means of circular annuli in the projection plane). Any emerging constraint on the morphology and degree of anisotropy of the stellar system should be, therefore, considered as resulting properties of the best-fit model, which has been selected exclusively on the basis of the (spherically averaged) kinematic and photometric information.

As for the morphological characterization (see Section 5.2), the projected isodensity contours are calculated on the basis of the nonspherical projected number density distribution. The relevant ellipticity profile is then constructed by considering the ratio of the principal axes of approximately 60 isodensity contours, corresponding to selected values of the normalized projected number density in the range $[0.9,{10}^{-3}];$ smooth profiles are then obtained by performing an average on subsets made of 10–20 individual ellipticity values.

For completeness, we also explored the distance value of 4.15 kpc from Watkins et al. (2015b, Paper III) in our fitting procedures. The best-fit model based on the 4.15 kpc value provides comparable results to those based on the 4.5 kpc distance value. However, the Paper III-based model offers a better fit to the LOS velocity dispersion profile, but a worse fit to the surface brightness profile.

5.1.1. Correction for Energy Equipartition Effects

Since PM measurements sample different kinematic tracers than LOS measurements—namely, stars with lower mass than bright red giant stars—some caution is needed when applying a one-component dynamical model simultaneously to the three-dimensional kinematics. In particular, our PM-based velocity dispersion profiles are constructed using MS stars within the magnitude range $20\lt {m}_{{\rm{F}}606{\rm{W}}}\lt 22$. This implies a typical stellar mass between 0.62 and 0.47 M_⊙ (using a Dotter et al. 2008 isochrone with [Fe/H] = −0.5, [α/Fe] = 0.2, d = 4.5 kpc, and E(B − V) = 0.04). The LOS measurements have instead a typical stellar mass of 0.83 M_⊙, since they sample only bright giant stars.

Since GCs reach a state of partial energy equipartition (see, e.g., Trenti & van der Marel 2013; Bianchini et al. 2016), their kinematics are expected to show a mass dependence, with lower-mass stars having gradually higher velocity dispersions. Bianchini et al. (2016) showed that the mass dependence of kinematics depends on the relaxation condition of the cluster (see their Equation (6)). Therefore, by knowing the relaxation state of a cluster it is possible to predict the shape of the velocity dispersion as a function of mass $\sigma (m)$ and then rescale the PM-based dispersion profiles according to the calculated factor.

Given ${n}_{\mathrm{rel}}={T}_{\mathrm{age}}/{T}_{{r}_{{\rm{c}}}}=169.6$ (number of core relaxation times the cluster has experienced), the mass dependence of the velocity dispersion from Equation (3) of Bianchini et al. (2016) is

$$\begin{eqnarray*}&&\sigma (m)={\sigma }_{0}\exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{m}{{m}_{\mathrm{eq}}}\right),\end{eqnarray*}$$

with ${m}_{\mathrm{eq}}=1.60$ M_⊙ (from Equation (6) of Bianchini et al. 2016) and ${\sigma }_{0}$ a normalization factor. Given a mass of 0.83 M_⊙ for the LOS velocities and a median value of 0.54 M_⊙ for the stars used for the PM-based velocity dispersion profiles (estimated via the adopted isochrone), the velocity dispersions of the two different mass tracers are related by σ(0.83)/σ(0.54) = 0.913. We rescale the tangential and radial PM dispersion profiles of the model by dividing them by this correction factor. We do not rescale the PM rotation profile, since we do not observe any signature of dependence of mass.

Note that the central velocity dispersion estimate based on the PM of evolved stars (Paper II) is 0.573 ± 0.005 $\mathrm{mas}\,{\mathrm{yr}}^{-1}$, which translates into 12.2 ± 0.1 $\mathrm{km}\,{{\rm{s}}}^{-1}$ at a distance of 4.5 kpc. The adopted rescaling factor of 0.913 implies a central velocity dispersion for RGB stars of 12.46 $\mathrm{km}\,{{\rm{s}}}^{-1}$, in full agreement with the value reported in Paper II.

5.2. Results

5.2.1. Projected Properties

The results of the fit are reported in Figures 6–9. The blue, green, and red lines in each plot correspond to the fit with inclination angles i = 25°, 30°, and 35°, respectively. The structural parameters of the three different models do not differ significantly; however, the model with i = 30° gives a better fit. We will consider this model as our fiducial best-fit model, and we report the best-fit parameters and structural properties in Tables 6 and 7.

Zoom In Zoom Out Reset image size

Table 6.  Dimensionless Parameters and Physical Scales of the Best-fit Model with $i=30^\circ$

Dimensionless Parameters			Physical Scales			Assumed Distance
Ψ	χ	$\bar{b}$	${{\rm{\Sigma }}}_{0}$	r0	v0	d
			(mag arcsec−2)	(arcsec)	(km s−1)	(kpc)
7.6	$1.6\times {10}^{-3}$	0.007	14.5 ± 0.1	26.3 ± 1.1	14.3 ± 0.1	4.5

Note. Concentration parameter Ψ, rotation strength parameter χ, the $\bar{b}$ parameter (the additional parameter c is set to unity), V-band central surface brightness ${{\rm{\Sigma }}}_{0}$ in mag arcsec⁻², radial scale r₀ in arcsec, velocity scale v₀ in km s⁻¹, and assumed distance d in kpc. For the physical scales, the associated $1\sigma$ errors are also shown.

Download table as:  ASCIITypeset image

Table 7.  Structural Parameters Derived for the Best-fit Model with i = 30°

C	${r}_{{\rm{c}}}$	${r}_{{\rm{h}}}$	${R}_{\mathrm{tr}}$	M	$M/{L}_{{\rm{V}}}$	$\mathrm{log}{\rho }_{0}$
	(arcsec)	(arcsec)	(arcsec)	(${M}_{\odot }$)	(M⊙/M⊙)	(M⊙ pc−3)
1.91 ± 0.02	26.5 ± 1.16	240.7 ± 10.6	1941.74 ± 85.67	$8.4\pm 0.4\times {10}^{5}$	1.98 ± 0.26	5.02 ± 0.01

Note. Concentration parameter $C=\mathrm{log}({R}_{\mathrm{tr}}/{r}_{{\rm{c}}})$, projected core radius ${r}_{{\rm{c}}}$ in arcsec, projected half-light radius ${r}_{{\rm{h}}}$ in arcsec, tridimensional truncation radius ${R}_{\mathrm{tr}}$ in cylindric coordinates (in arcsec), total mass of the cluster M in units of M_⊙, V-band mass-to-light ratio in solar units, logarithm of the central mass density ${\rho }_{0}$ in units of M_⊙ pc⁻³.

Download table as:  ASCIITypeset image

Figure 6 shows the rotation profile based in the LOS direction (right) and in the plane of the sky (left). The additional yellow lines represent the best-fit model derived by Bianchini et al. (2013) assuming an inclination of 45° and without accounting for the rotation in the plane of the sky. The figure clearly shows that our new best-fit model, assuming the new inclination value of i ≃ 30°, is able to reproduce the three-dimensional rotational structure of 47 Tuc, whereas a higher inclination angle would fail in describing the rotation in the plane of the sky.

Our model predicts an anisotropy profile that is in excellent agreement with the observations (Figure 7), characterized by isotropy in the central cluster regions and mild radial anisotropy in the intermediate regions. Moreover, we are able to simultaneously reproduce all the dispersion profiles (Figure 8) and the surface brightness profile (Figure 9), employed for the fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the adopted dynamical equilibria are axisymmetric, we can calculate the corresponding projected ellipticity profile, which we illustrate in Figure 10, together with the ellipticity data currently available for 47 Tuc (from White & Shawl 1987). We recall that the ellipticity profile associated with the best-fit self-consistent model is a structural property completely determined by the values of the dimensionless parameters and physical scales identified during the model selection procedure. As already appreciated in our previous analysis (see Bianchini et al. 2013), we find a very good agreement with the available observational profile. Such a result allows us to confirm, with increased confidence, that the physical origin of the observed flattening of 47 Tuc is indeed the presence of internal rotation. We emphasize that such an agreement is nontrivial, especially in consideration of the reduced value of the inclination angle (in our previous study we adopted i = 45°, while now we determined i = 30° to be more appropriate). The appreciable morphological consistency recovered, once again, in the current analysis should be interpreted as a manifestation of the more significant intrinsic rotation (as recovered from the additional constraints posed by new PM data sets), which compensates the effects of a less favorable LOS direction. For completeness, in Figure 11 we also show the projected isodensity contours of the discrete realization of the best-fit axisymmetric model as a function of the tridimensional radius in cylindrical coordinates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2.2. Intrinsic Properties

Our three-dimensional model allows us to explore the intrinsic kinematic and morphological structure of 47 Tuc. The derived intrinsic V/σ profile, measured in the equatorial plane perpendicular to the rotation axis, is characterized by a peak value of ∼0.9 reached at around two half-light radii (Figure 12). This confirms that the internal rotation of 47 Tuc is higher than what has been reported in previous studies (e.g., the LOS-based value V/σ = 0.25; Bianchini et al. 2013). Correspondingly, we can also characterize the three-dimensional structure of the cluster by means of the intrinsic ellipticity profile, as measured in the meridional plane of the system (defined by the rotation axis and any of the principal axes on the equatorial plane) and expressed as a function of the semimajor axis (Figure 13). The behavior of this profile in the central regions is approximately linear, as shaped by the solid-body-like behavior of the rotation curve in that portion of the system. We emphasize that the radial location of the peak of the ellipticity profile does not correspond to the location of the peak of the intrinsic V/σ profile, as appropriate in the presence of a nontrivial coupling between the angular momentum and mass distribution within the system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We recall that the intrinsic morphological and kinematic structure of an axisymmetric (nonstratified) rotating equilibrium is more complex than a simple characterization expressed by means of unidimensional radial profiles; therefore, we include also a representation of the bidimensional maps depicting the structure of the density and angular momentum distribution, expressed in terms of the contours of the normalized isodensity and isovelocity surfaces evaluated on the meridional plane of the three-dimensional model. In the maps illustrated in Figure 14 the radial range on the equatorial plane is approximately equal to $2{r}_{{\rm{h}}}$. From the density map (left panel) it may be appreciated that the maximum degree of flattening is reached in the intermediate regions of the system and that the shape of the isodensity contours becomes more complex as the density decreases (i.e., the outer contours are not monotonically increasing functions of the polar angle, hence the dimples on the rotation axis). From the kinematic map (right panel) it may be noted that the velocity field is not cylindrically stratified, presenting an absolute maximum at about 300''; the location of the peak of the intrinsic V/σ profile is determined by the interplay between the three-dimensional structure of such a velocity field and that of the trace of the velocity dispersion tensor. Equivalent maps may be calculated also for the meridional sections of the equipotential and isobaric surfaces.

Zoom In Zoom Out Reset image size

Lastly, we wish to emphasize that the current presence of such an appreciable degree of internal rotation in 47 Tuc carries important implications about the initial amount of angular momentum of the cluster. Indeed, many facets of the role played by internal rotation during the long-term evolution of collisional stellar systems must still be explored, but there is confirmed evidence that two-body relaxation and mass loss determine transport and loss of angular momentum in clusters (see, e.g., Einsel & Spurzem 1999; Ernst et al. 2007); therefore, any measurement of the presence of internal rotation at the present day should be considered as a lower limit of the initial angular momentum content. The emerging kinematic complexity of 47 Tuc (and a progressively increasing number of Galactic GCs) therefore offers an essential ingredient toward a more complete and fundamental understanding of the formation and early dynamical evolution of GCs.

5.2.3. Some Remarks on the Modeling Strategy

We recognize that our modeling strategy has a number of limitations, especially regarding the description of the projected structural and kinematic profiles in the outer regions of the cluster. With particular reference to the behavior of the surface brightness profile and LOS velocity dispersion profile in the very outer parts (at radii >1500''), we interpret their discrepancies as due to the effects of the tidal field of the Milky Way on the structure and internal kinematics of the system. It is well known that, during the course of their dynamical evolution (as driven by two-body relaxation), collisional stellar systems tend to expand until they fill their Roche lobe, and then they progressively start to lose mass. Such a process significantly affects the phase-space structure of their periphery, both by increasing the complexity of the kinematics of the stars that are energetically bound (e.g., Tiongco et al. 2016a, 2017) and by determining the existence of a population of energetically unbound stars that are nonetheless spatially confined within the cluster ("potential escapers"; see Küpper et al. 2010; Claydon et al. 2017). In terms of projected observables, these evolutionary effects may manifest themselves in the form of surface brightness profiles extending beyond the cutoff radius predicted by a simple spherical King model (i.e., the so-called "extra-tidal" structures identified in many globular clusters in Local Group galaxies; see, e.g., Barmby et al. 2002; Harris et al. 2002; McLaughlin & van der Marel 2005) and velocity dispersion profiles characterized by an untruncated, flattened behavior (see, e.g., the studies of M15 and M92, Drukier et al. 1998, 2007; and NGC 5694, Bellazzini et al. 2015). In term of intrinsic properties, tidally perturbed systems are also characterized by a flavor and degree of anisotropy that strongly depends on the strength of the tidal environment in which they have evolved (see Tiongco et al. 2016b, and additional remarks at the end of next paragraph).

In this respect, we emphasize that simple, physically based dynamical models of the kind adopted in our study, as defined by a quasi-Maxwellian distribution function, suitably modified near the tidal boundary and truncated above it, successfully capture the phase-space properties of the bulk of the cluster members by relying on a truncation prescription that heuristically mimics the effects of the tidal field (e.g., Woolley 1954; King 1966; Wilson 1975; and see, more recently, Gieles & Zocchi 2015). Of course, the choice of such a one-dimensional (energy) truncation prescription in the definition of the distribution function strongly affects the structural and kinematic properties of the resulting configurations (see Davoust 1977; Hunter 1977), but, by definition, they cannot offer a realistic description of the tidal field (see Heggie & Ramamani 1995; Bertin & Varri 2008) or account for the existence of energetically unbound members (see Daniel et al. 2017). We wish to emphasize that, despite their inherent simplicity, the self-consistent rotating models adopted here have nonetheless two advantageous properties. First, their phase-space truncation reduces, in the nonrotating limit, to the smooth Wilson (1975) prescription, which is now often considered more successful in describing the outer regions of clusters than the traditional King (1966) model (see McLaughlin & van der Marel 2005), which demands only continuity in phase space. Second, their velocity dispersion tensor is characterized by isotropy in the central region, weak radial anisotropy in the intermediate regions, and tangential anisotropy in the outer parts. Such a behavior of the pressure tensor was not assigned a priori in the definition of the models, but it results from the requirement of self-consistency, positivity of the distribution function, and energy truncation in the presence of angular momentum as a second integral of the motion. Isotropy or mild tangentially biased pressure anisotropy in the outer parts of a star cluster has been shown to be a natural result of the dynamical evolution of a collisional stellar system within an external tidal field, which induces a preferential loss of stars on radial orbits (see, e.g., Giersz & Heggie 1997; Takahashi et al. 1997; Baumgardt & Makino 2003; Vesperini et al. 2014; Tiongco et al. 2016b; Claydon et al. 2017).

The other discrepancy, namely, a model underestimate of the surface brightness in the very central region and in an intermediate region in the range 100''–200'', can instead be considered as a result of a χ²-minimization employed by the fitting procedure that gives more weight to the three-dimensional kinematic data (see Section 5.1). As a result, the value of the best-fit half-light radius we obtain (${r}_{h}=240.7\pm 10\buildrel{\prime\prime}\over{.} 9$) is larger than what is reported in previous literature works (e.g., ${r}_{h}=190\buildrel{\prime\prime}\over{.} 2$; Harris 1996).

We wish to emphasize that our study focuses on understanding the rotational structure of 47 Tuc. The kinematic profiles (the two rotation curves and the three velocity dispersion profiles) show a general agreement within the observational error bars. However, our model underestimates the values of the rotation profiles in the intermediate region and overestimates the value of velocity dispersion along the LOS. These discrepancies could be interpreted in view of the intrinsic limitations of the choice of the distribution function defining the dynamical models under consideration. Indeed, alternative modeling strategies may offer additional degrees of freedom in the parameterization of the phase space, but we wish to note that differentially rotating equilibria are relatively rare, especially for globular cluster studies (see Sollima et al. (2009) for an application of the rotating models proposed by Wilson (1975) to the study of ω Cen). An alternative approach for the construction of phase-space equilibria is based on the use of actions instead of integrals of the motion; an example of this line of attack has been recently proposed by Jeffreson et al. (2017). Finally, we stress that the velocity anisotropy profile and the ellipticity profile are not directly used in the fitting procedure, and the respective figures show a comparison between the available data and the model prediction. This result may also be interpreted as a reassuring "a posteriori" validation of our choice of adopting axial symmetry and reduces the scope of exploring more general symmetries (e.g., triaxial configurations).

A complementary path would consist in employing empirical models that are optimally designed to describe the data (e.g., Jeans models or orbit-based Schwarzschild models). Despite providing more freedom in the description of the data, the major drawback is that they do offer a very limited connection to the underlying physical picture of the stellar system under consideration. However, for studies aimed at the understanding of specific physical ingredients (e.g., the possible presence of an intermediate-mass black hole in 47 Tuc, the exact mass-to-light ratio of the stellar population, or the amount of energy equipartition), we would certainly benefit from complementary and alternative modeling approaches from the one employed here. In addition, thanks to recent progress on the computational side, it has recently been proved that realistic N-body models, following the entire dynamical evolution of the system, are finally within reach, at least for selected Galactic globular clusters (see, e.g., the model of M4 presented by Heggie 2014 or the DRAGON series by Wang et al. 2015), although it should be kept in mind that the performance cost of the exploration of a wide range of initial conditions is still far from negligible, especially for massive clusters such as 47 Tuc.

In order to understand the limitations intrinsic to our modeling, we report in Table 8 a compilation of the total mass and the mass-to-light ratio for 47 Tuc from the literature obtained using a variety of dynamical modeling techniques. These include (i) distribution-function-based models (spherical and isotropic models, McLaughlin & van der Marel 2005; Zocchi et al. 2012; spherical and anisotropic models, de Vita et al. 2016; axisymmetric rotating models, Bianchini et al. 2013), (ii) isotropic spherical Jeans models (Watkins et al. 2015a), and (iii) N-body and spherical Monte Carlo models (Giersz & Heggie 2011; Baumgardt 2017). Note that all these previous works, with the exception of Bianchini et al. (2013), do not take into consideration internal rotation; therefore, they are intrinsically unable to reproduce our state-of-the-art three-dimensional kinematic measurements.

Table 8.  Comparison with Previous Works

M	$M/{L}_{{\rm{V}}}$	Dynamical Model	References
(105 M⊙)	(M⊙/M⊙)
10.7 ± 0.98	${1.17}_{-0.43}^{+0.53}$	Spherical Wilson (1975)	McLaughlin & van der Marel (2005)
9.0	1.52	Spherical Monte Carlo model	Giersz & Heggie (2011)
7.18 ± 0.41	1.34 ± 0.08	Spherical King (1966)	Zocchi et al. (2012)
6.23 ± 0.04	1.69 ± 0.13	Rotating Varri & Bertin (2012)	Bianchini et al. (2013)
${5.57}_{-0.28}^{+0.33}$	1.40 ± 0.03	Isotropic Jeans model	Watkins et al. (2015a)
7.00 ± 0.06	1.99 ± 0.20	N-body model	Baumgardt (2017)
7.77	1.63	Anisotropic ${f}_{T}^{(\nu )}$ model	de Vita et al. (2016)
8.4 ± 0.4	1.98 ± 0.26	Rotating Varri & Bertin (2012)	This work

Note. Compilation of the total mass and the mass-to-light ratio of 47 Tuc derived from dynamical modeling in previous works.

Download table as:  ASCIITypeset image

From Table 8, it is evident that the estimates in the literature depend on the particular modeling techniques employed. The values of the total mass of 47 Tuc are in the range of (5.5–10.7) × 10⁵ M_⊙, with our model giving an intermediate value of 8.4 × 10⁵ M_⊙, while the mass-to-light ratio is in the range of 1.2–2 M_⊙/ L_⊙, with our model giving 1.98 M_⊙/ L_⊙, consistent with the upper limit. Given the agreement with these previous works, we are confident that our model is able to reproduce the global properties of 47 Tuc, despite some of the discrepancies reported above.

We have derived for the first time the plane-of-the-sky rotation of the GC 47 Tuc from the core out to 13' (or three to four half-light radii, or ∼30% of the tidal radius), together with detailed radial and tangential velocity dispersion profiles. The cluster's plane-of-the-sky kinematic information is coupled to literature LOS measurements of the same quantities and surface brightness profiles and simultaneously fit with state-of-the-art dynamical models. From the application of the rotating dynamical model we conclude the following:

• 1.  
  PMs are critically necessary if we hope to constrain the intrinsic structure of GCs. With the application of our dynamical model to both kinematics and photometry, we have obtained a full three-dimensional description of 47 Tuc.
• 2.  
  The higher rotation in the plane of sky with respect to the one along the LOS implies an inclination angle of the rotation axis of ≃30°.
• 3.  
  Our best-fit dynamical model predicts an intrinsic V/σ profile that reaches values of ≃0.9 around two half-light radii from the cluster's center.
• 4.  
  On the basis of our global dynamical analysis, we confirm that the observed flattening of the cluster is most likely due to its appreciable internal rotation, as the projected ellipticity profile determined by our model is in good agreement with the ellipticity measurements currently available for 47 Tuc. The three-dimensional morphological structure implied by our axisymmetric model is complex and may be characterized by a nonmonotonic intrinsic ellipticity profile, which reaches values of about 0.45 toward the intermediate regions of the cluster.

This comprehensive dynamical investigation, based on new PM measurements of unprecedented accuracy, has allowed us to unveil a new degree of kinematic complexity in 47 Tuc. Such a superb characterization of the three-dimensional velocity space offers a novel and invaluable ground for the study of numerous aspects of collisional gravitational dynamics, from the long-forgotten role played by angular momentum to the tantalizing opportunity of finally exploring the phase-space properties of its stellar populations, with the goal of providing a key contribution toward a more realistic dynamical paradigm for this class of stellar systems.

A.B. acknowledges support from STScI grant AR-12845 (PI: Bellini). P.B. acknowledges financial support from CITA National Fellowship. A.L.V. acknowledges support from the EU Horizon 2020 program in the form of a Marie Sklodowska-Curie Fellowship (MSCA-IF-EF-RI 658088). G.P. acknowledges partial support by PRIN-INAF 2014 and by the "Progetto di Ateneo 2014 CPDA141214" by Università di Padova.