Hubble Space Telescope Proper Motion (HSTPROMO) Catalogs of Galactic Globular Clusters. V. The rapid rotation of 47 Tuc traced and modeled in three dimensions

TL;DR

This study uses high-precision proper motions to model the three-dimensional rotation and anisotropy of globular cluster 47 Tuc, revealing complex internal dynamics and emphasizing the importance of proper motions in understanding cluster structure.

Contribution

It introduces a comprehensive 3D dynamical model of 47 Tuc incorporating proper motions, rotation, and anisotropy, improving upon previous line-of-sight only analyses.

Findings

Measured the cluster's rotation in the plane of the sky for the first time.

Revealed a peak V/sigma ratio of ~0.9 at two half-light radii.

Discovered a non-monotonic ellipticity profile up to 0.45.

Abstract

High-precision proper motions of the globular cluster 47 Tuc have allowed us to measure for the first time the cluster rotation in the plane of the sky and the velocity anisotropy profile from the cluster core out to about 13'. These profiles are coupled with prior measurements along the line of sight and the surface-brightness profile, and fit all together with self-consistent models specifically constructed to describe quasi-relaxed stellar systems with realistic differential rotation, axisymmetry and pressure anisotropy. The best-fit model provides an inclination angle i between the rotation axis and the line-of-sight direction of 30 deg, and is able to simultaneously reproduce the full three-dimensional kinematics and structure of the cluster, while preserving a good agreement with the projected morphology. Literature models based solely on line-of-sight measurements imply a significantly different inclination angle (i=45 deg), demonstrating that proper motions play a key role in constraining the intrinsic structure of 47 Tuc. Our best-fit global dynamical model implies an internal rotation higher than previous studies have shown, and suggests a peak of the intrinsic V/sigma ratio of ~0.9 at around two half-light radii, with a non-monotonic intrinsic ellipticity profile reaching values up to 0.45. Our study unveils a new degree of dynamical complexity in 47 Tuc, which may be leveraged to provide new insights into the formation and evolution of globular clusters.

Tables

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\times 765s$, $1\times 1200s$
			F814W	$1\times 690s$, $1\times 1020s$
2006.67	10730	ACS/WFC	F435W	$8\times 350s$
	10737	ACS/WFC	F435W	$1\times 339s$
			F555W	$1\times 339s$
			F606W	$1\times 339s$
			F814W	$1\times 339s$
2009.54	11444	WFC3/UVIS	F606W	$24\times 350s$
	11452	WFC3/UVIS	F814W	$6\times 350s$
2012.22	12692	WFC3/UVIS	F606W	$2\times 350s$
2014.35	13596	ACS/WFC	F435W	$2\times 339s$
			F475W	$2\times 339s$
			F606W	$2\times 339s$
			F775W	$2\times 339s$
			F814W	$2\times 339s$
			F850L	$2\times 339s$

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\times 600s$, $7\times 1000s$
			F814W	$7\times 1000s$
1994.86	5369	WFPC2	F606W	$9\times 300s$, $1\times 500s$
1999.79	8095	WFPC2	F606W	$1\times 400s$, $7\times 500s$
2013.34	12971	ACS/WFC	F435W	$1\times 290s$, $1\times 690s$
			F555W	$1\times 360s$, $1\times 660s$

Table 3
List of Observations of the Outer Field
Epoch	GO	Instr.	Filter	$N\times \exp .\mathrm{time}$
2002.35	9318	WFPC2	F606W	$1\times 300s$, $1\times 400s$
				$7\times 500s$
2002.40	8059	WFPC2	F450W	$1\times 300s$
			F606W	$1\times 160s$
			F814W	$1\times 100s$,$1\times 300s$
2003.40	9709	WFPC2	F606W	$1\times 400s$,$7\times 500s$
2010.34	11677	WFC3/UVIS	F390W	$1\times 1048s$,$1\times 1099s$
				$1\times 1212s$,$1\times 1355s$
				$2\times 1400s$
			F606W	$1\times 1252s$,$3\times 1347s$
				$1\times 1398s$,$1\times 1402s$

Table 4
Rotation of 47 Tuc in the Plane of the Sky
$⟨r⟩$	NSTAR	$⟨{\mu }_{\mathrm{TAN}}⟩$	$⟨{\sigma }_{{\mu }_{\mathrm{TAN}}}⟩$	$⟨r⟩$	NSTAR	$⟨{\mu }_{\mathrm{TAN}}⟩$	$⟨{\sigma }_{{\mu }_{\mathrm{TAN}}}⟩$
(′′)		($masyr$)	($\mathrm{mas}{\mathrm{yr}}^{-1}$)	(′′)		($\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)†		($-$0.233)†	(0.055)†	705.53	118	$-$0.232	0.041
352.61	150	$-$0.283	0.016	775.05	117	$-$0.208	0.045
$†$ From Anderson & King (2003a).

Table 5
Radial and tangential velocity-dispersion
profile data of 47 Tuc
$⟨r⟩$	NSTAR	$⟨{\sigma }_{\mathrm{RAD}}⟩$	err${}_{{\sigma }_{\mathrm{RAD}}}$	$⟨{\sigma }_{\mathrm{TAN}}⟩$	err${}_{{\sigma }_{\mathrm{TAN}}}$
(′′)		($\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

Table 6
Dimensionless Parameters and Physical Scales of the Best-Fit Model with $i={30}^{\circ }$
	Dimensionless parameters				Physical scales				Assumed distance
	$\Psi$	$\chi$	$\overline{b}$		${\Sigma }_0$	$r_0$	$v_0$		$d$
					mag arcsec-2	arcsec	km s-1		kpc
	7.6	$1.6\times {10}^{-3}$	0.007		$14.5\pm 0.1$	$26.3\pm 1.1$	$14.3\pm 0.1$		$4.5$

Table 7
Structural Parameters Derived for the Best-Fit Model with $i={30}^{\circ }$
$C$	$r_c$	$r_h$	$R_{\mathrm{tr}}$	$M$	$M/L_V$	$\log {\rho }_0$
	arcsec	arcsec	arcsec	$M_{\odot }$	$M_{\odot }/L_{\odot }$	$M_{\odot }$  pc-3
$1.91\pm 0.02$	$26.5\pm 1.16$	$240.7\pm 10.6$	$1941.74\pm 85.67$	$8.4\pm 0.4\times {10}^5$	$1.98\pm 0.26$	$5.02\pm 0.01$

Table 8
Comparison with Previous Works
$M$	$M/L_V$	dynamical model	reference
(${10}^5$ $M_{\odot }$)	($M_{\odot }/L_{\odot }$)
$10.7\pm 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\pm 0.41$	$1.34\pm 0.08$	spherical King (1966)	Zocchi et al. (2012)
$6.23\pm 0.04$	$1.69\pm 0.13$	rotating Varri & Bertin (2012)	Bianchini et al. (2013)
${5.57}_{-0.28}^{+0.33}$	$1.40\pm 0.03$	isotropic Jeans model	Watkins et al. (2015a)
$7.00\pm 0.06$	$1.99\pm 0.20$	$N$-body model	Baumgardt (2017)
$7.77$	$1.63$	anisotropic $f_T^{\left(\nu \right)}$ model	de Vita et al. (2016)
$8.4\pm 0.4$	$1.98\pm 0.26$	rotating Varri & Bertin (2012)	this work

Equations

$$\sigma (m)={\sigma }_0\exp \left(-\frac{1}{2}\frac{m}{m_{\mathrm{eq}}}\right),$$

Full text

1. Introduction

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 to calculate 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; van Leeuwen et al. 2000; Meylan & Heggie 1997 and references therein), 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., Noyola et al. 2010; Anderson & van der Marel 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 with 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 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. 2016a, 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. 2016b).

On the observational, proper-motion-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 $\omega$ Centauri (see their Fig. 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 $\sim$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 4 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 the following: in Sect. 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 Sect. 3 and 4. In Sect. 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 Sect.6.

2. Data sets and Reduction

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 arcsecs (>$$20^{\prime}) on the sky. In Fig. 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 the East to the West, as: the inner field, the central field, the calibration field, and the outer field. We made use of exposures taken through 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).

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 PSFs (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); Bellini, Anderson & Bedin (2011). The photometric data of the central field comes 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 six 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, 6-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

The PM catalog of the central field is that published in Paper I. We analyzed 433 exposures taken with the ACS (both HRC and WFC channels) and the WFC3/UVIS. The complete list of observations can be found in Table 7 of Paper I. The available time baseline used to compute the motion of each star goes from 0.52 to 10.32 years (median of 8.23), depending on the available overlap between different exposures. We have 103 638 stars with measured PMs in the central field. We closely applied the prescriptions given in Section 7.5 of Paper I in order to select only high-quality PM measurements (53 898 stars).

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.

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 Declination direction, and all the available WFPC2 exposures overlapping them. The complete list of observations for the inner and outer fields are reported in Tables 2 and 3, respectively.

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 that what we can achieve with later HST optical imagers. The WFPC2 camera was composed of four 800$\times$800 pixel2 12-bit detectors, with a pixel scale of 99.6 mas pixel*-1* for the three WF chips, and 45.5 mas pixel*-1* for the PC chip. All four chips suffered from $\sim$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 years, 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 to minimize inter-chip 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 Sect. 5.2 of Paper I for more details).

The computed PMs for the inner and the outer fields are in agreement with with those computed in the central and calibration fields in terms of expected intrinsic values and associated errors. On the other hand, the available dataset 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 the subtraction in quadrature of PM errors) using inner- and outer-field measurements as for the central and calibration field.

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.

3. The rotation of 47 Tuc in the plane of the sky

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.

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 Fig. 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 counter-clockwise 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 Fig. 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 sigma outside the distribution of 47 Tuc stars. This rough selection necessarily includes field stars and a few cluster members. We then iteratively computed 4$\sigma$-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 $\sigma$ within 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 kept only SMC stars within 4 $\sigma$ from their median location.

The final sample of SMC stars is what we show in red in Fig. 2. We have, respectively: 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_{\rm 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$\sigma$-clipped estimate of their barycenter location on the PM diagram. This location has coordinates: $\mu_{\alpha}\cos\delta=-4.797$ ${\rm mas}\,{\rm yr}^{-1}$, $\mu_{\delta}=1.244$ ${\rm mas}\,{\rm yr}^{-1}$, and it defines the center of the crosshairs in the middle panels of Fig. 2.

We zero-pointed the PM of SMC stars to the reference point defined above, and computed their average $\mu_{\rm 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 2 measurements for the outer field (117 stars each). The computed quantities, as a function of the radial distance, are shown with errorbars in Fig. 3, and are listed in Table 4. We adopted the convention of using negative $\mu_{\rm TAN}$ values in case of a clockwise rotation in the plane of the sky (North to West). The horizontal bars in Fig. 3 indicate the size of the radial intervals within which each $\mu_{\rm 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_{\rm TAN}=mr$ gives $m=-1.65\pm 0.10$ ${\rm mas}\,{\rm yr}^{-1}$ arcsec*-1*. Then, the rotation profile slowly flatters at larger radii, to reach the highest value of $-0.312\pm 0.015$ ${\rm mas}\,{\rm yr}^{-1}$ at $\sim 390^{\prime\prime}$ ($\sim 6\farcm 52$) from the cluster center. The measured rotation slows down to about $-0.21$ ${\rm mas}\,{\rm yr}^{-1}$ in the outermost regions probed by our data.

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 Sect. 5.

To give the reader a better sense of how the cluster is rotating on the plane of the sky, we show in Fig. 4 the rotation map of the cluster derived with the datasets we analyzed. To obtain the map, we divided the datasets 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 ${\rm mas}\,{\rm yr}^{-1}$, has been magnified by a factor 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).

4. Velocity-dispersion and anisotropy profiles

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 field, to become again somewhat more circular in the outer field (middle panels of Fig. 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 datasets 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<$$m_{\rm F606W}$$<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 Fig. 5 shows the radial velocity-dispersion profile $\sigma_{\rm RAD}$ of 47 Tuc as a function of radius. In red we report the values computed in the inner and outer field, for completeness. The horizontal bars illustrate the radial extent over which each point is obtained. The central radial velocity dispersion is about 0.64 ${\rm mas}\,{\rm yr}^{-1}$, or 13.65 ${\rm km}\,{\rm s}^{-1}$ at a distance of 4.5 kpc (to be compared to 11.5 ${\rm km}\,{\rm s}^{-1}$ for evolved stars, Harris 1996).111We are assuming here that the LOS-based velocity-dispersion-profile value quoted in the Harris catalog, which is the average of available measurements in the literature, refers to the center of the cluster. This might actually not be true, as literature values are computed over different radial ranges. We will see later, in Section 5.1.1, that our $\sigma$-mass dependence modeling predicts a red-giant-branch- (RGB-) mass scaled central velocity-dispersion value closer to about 12.5 ${\rm km}\,{\rm s}^{-1}$ rather than 11.5 ${\rm km}\,{\rm s}^{-1}$. The profile reaches a minimum value of about 0.31 ${\rm mas}\,{\rm yr}^{-1}$ in the outermost point of the calibration-field (\sim$$8\farcm 25, or about 495*′′*). The middle panel of the figure shows the tangential velocity-dispersion profile $\sigma_{\rm TAN}$ computed in the same radial intervals as for $\sigma_{\rm RAD}$. The two vertical lines mark the location of the core radius $r_{\rm c}$=$0\farcm 36$ ($21\farcs 6$) and the half-light radius $r_{\rm h}$=$3\farcm 17$ ($190\farcs 2$) (Harris 1996). Table 5 lists the values of the data points shown in Fig. 5.

The deviation from isotropy ($\sigma_{\rm TAN}$/$\sigma_{\rm RAD}$-1) as a function of radius is shown in the bottom panel of Fig. 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 the second points of the dispersion profiles are nonetheless consistent with each other to within less than one sigma. Changing the size of the first few radial bins produced similar trends, all consistent with being flat in the centermost regions within the errorbars. Also note that these profiles are based on relatively faint 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 S/N RGB stars, does not show any central drop.

5. Dynamical models

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 (Kacharov et al., 2014, Bianchini et al., 2013).

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 $\Psi$, rotation strength parameter $\chi$, 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 single-mass component, they do not take into consideration the effects connected with mass segregation. The implications for this assumption are described in Sect. 5.1.1.

We will use the PM-based profiles described in Sect. 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 Sect. 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_{0}$, the central surface density $\Sigma_{0}$, and the velocity scale $v_{0}$), 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^{\circ}$, derived by applying Eqn. 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^{\circ}$–$35^{\circ}$). Since PMs are measured in mas yr*-1* in the plane of the sky, a multiplying factor of $4.74\,d$ is needed to convert these values to ${\rm km}\,{\rm s}^{-1}$. We will assume the distance of 47 Tuc of $d=4.5$ kpc (Harris 1996).

The fitting procedure is twofold. First, we determine the dimensionless parameters needed to reproduce the observed value of $V/\sigma$ and the observed position of the rotation peak (for further details see Sect. 3.1 and 3.5 of Bianchini et al. 2013). Then we calculate the physical scales through $\chi^{2}$ minimization for all the kinematic profiles (three dispersion profiles, two rotation curves) to obtain the radial scale $r_{0}$ and the velocity scale $v_{0}$, and then we fit the surface-density profile to obtain the central surface density $\Sigma_{0}$. This provides all the constrains 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 by 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 to the morphological characterization (see Sect. 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<$$m_{\rm F606W}$$<22. This implies a typical stellar mass between 0.62 and 0.47 $M_{\odot}$ (using a Dotter et al. 2008 isochrone with [Fe/H]=$-$0.5, [$\alpha$/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_{\odot}$, 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 Eqn. 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_{\rm rel}=T_{\rm age}/T_{r_{\rm c}}=169.6$ (number of core relaxation times the cluster has experienced), the mass dependence of the velocity dispersion from Eqn. 3 of Bianchini et al. (2016) is

[TABLE]

with $m_{\rm eq}=1.60$ $M_{\odot}$ (from Eqn. 6 of Bianchini et al., 2016) and $\sigma_{0}$ a normalization factor. Given a mass of 0.83 $M_{\odot}$ for the LOS velocities and a median value of 0.54 $M_{\odot}$ 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 $\sigma(0.83)/\sigma(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$\pm$0.005 ${\rm mas}\,{\rm yr}^{-1}$, which translates into 12.2$\pm$0.1 ${\rm 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 ${\rm 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 Figs. 6–9. The azure, 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^{\circ}$ 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.

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\simeq 30^{\circ}$, 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 (Fig. 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 (Fig. 8) and the surface-brightness profile (Fig. 9), employed for the fit.

Since the adopted dynamical equilibria are axisymmetric, we can calculate the corresponding projected ellipticity profile, which we illustrate in Fig. 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^{\circ}$, while now we determined $i=30^{\circ}$ 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 datasets), which compensates the effects of a less favorable LOS direction. For completeness, in Fig. 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.

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/$\sigma$ profile, measured in the equatorial plane perpendicular to the rotation axis, is characterized by a peak value of $\sim$0.9 reached at around 2 half-light radii (Fig. 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/$\sigma$=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 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 (Fig. 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/\sigma$ profile, as appropriate in the presence of a nontrivial coupling between the angular momentum and mass distribution within the system.

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 Fig. 14 the radial range on the equatorial plane is approximately equal to $2r_{\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/\sigma$ 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.

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 towards 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 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. 2016b, 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 cut-off 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. McLaughlin & van der Marel 2005; Barmby et al. 2002; Harris et al. 2002) and velocity dispersion profiles characterized by an untruncated, flattened behavior (e.g., see the studies of M15 and M92 by Drukier et al. 1998, 2007 and NGC 5694 by 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. 2016a 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 mimic the effects of the tidal field (e.g. Woolley 1954; King 1966; Wilson 1975 and, 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 can not 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 non-rotating 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 (e.g., see Takahashi et al. 1997; Giersz & Heggie 1997; Baumgardt & Makino 2003; Claydon et al. 2017; Tiongco et al. 2016a; Vesperini et al. 2014).

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 $\chi^{2}$-minimization employed by the fitting procedure that gives more weight to the three-dimensional kinematic data (see Sect. 5.1). As a result, the value of the best-fit half-light radius we obtain ($r_{h}=240.7\pm 10.9$ arcsec) is larger than what is reported in previous literature works (e.g., $r_{h}=190.2$ arcsec,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 parametrization 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 $\omega$ 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 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 a 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 (e.g., see 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, and 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 (Baumgardt 2017; Giersz & Heggie 2011). 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.

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 ranges between 5.5–10.7$\times$105 $M_{\odot}$ with our model giving an intermediate value of 8.4$\times$105 $M_{\odot}$, while the mass-to-light ratio ranges from 1.2–2 $M_{\odot}/L_{\odot}$, with our model giving 1.98 $M_{\odot}/L_{\odot}$, 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.

6. Conclusions

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 3–4 half-light radii, or $\sim$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 that:

• •

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;

• •

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 $\simeq 30^{\circ}$;

• •

Our best-fit dynamical model predicts an intrinsic V/$\sigma$ profile that reaches values of $\simeq 0.9$ around two half-light radii from the cluster’s center;

• •

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 non-monotonic intrinsic ellipticity profile, which reaches values of about 0.45 towards 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 aspect 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 towards a more realistic dynamical paradigm for this class of stellar systems

AB acknowledges support from STScI grant AR-12845 (PI: Bellini). PB acknowledges financial support from CITA National Fellowship. ALV acknowledges support from the EU Horizon 2020 program in the form of a Marie Sklodowska-Curie Fellowship (MSCA-IF-EF-RI 658088). GP acknowledges partial support by PRIN-INAF 2014 and by the “Progetto di Ateneo 2014 CPDA141214” by Università di Padova.