© ESO, 2012

1. Introduction

A large class of stellar systems is expected to be in a quasi-equilibrium state characterized by a distribution function not far from a Maxwellian. Mechanisms that are thought to operate in this direction range from standard two-body gravitational scattering (Chandrasekhar 1943) to violent relaxation, instabilities, and phase mixing for collisionless systems (Lynden-Bell 1967). In particular, for relatively small stellar systems, such as globular clusters, two-star collisional relaxation is thought to be responsible for such quasi-relaxed condition, while for large stellar systems, such as elliptical galaxies, incomplete violent relaxation has likely acted during the formation process so as to bring the system to a state in which significant (radially biased) anisotropy is present only for weakly bound stars. To be sure, full relaxation would lead to a singular state, characterized by infinite mass; the so-called isothermal sphere solution may provide an approximate representation only of the inner parts of real stellar systems.

In line with the above considerations, several physically motivated self-consistent finite-mass models (thus characterized by significant, but not arbitrary, deviations from a pure Maxwellian) have been constructed and studied, with interesting astronomical applications. For the description of partially relaxed systems, self-consistent anisotropic models have been introduced to describe the products of incomplete violent relaxation (see Bertin & Stiavelli 1993; Trenti et al. 2005, and references therein). For the description of quasi-relaxed globular clusters, the spherical isotropic King (1966) models, which incorporate in a heuristic way the concept of tidal truncation, provide a standard reference framework (see Zocchi et al. 2012, for a photometric and kinematic study of a sample of Galactic globular clusters under different relaxation conditions). Recently, King models have been extended to the full triaxial case, to include explicitly the three-dimensional role of an external tidal field in the simple case of a circular orbit of the cluster within the host galaxy (Bertin & Varri 2008; Varri & Bertin 2009, in the following denoted as Papers I and II, respectively; see also Heggie & Ramamani 1995); these triaxial models can thus be seen as a collisionless analogues of the purely tidal Roche ellipsoids.

For astronomical applications, one of the main drivers (or empirical clues) in the above modeling investigations is the issue of the origin of the observed geometry of stellar systems, being it known that pressure anisotropy, tides, or rotation can be responsible, separately, for deviations from spherical symmetry. In the above-mentioned studies, very little attention was actually placed on the role of rotation. For ellipticals, most of the attention that led to the development of stellar dynamical models, after the first kinematical measurements became available in the mid-70s, was taken by the study of the curious behavior of pressure-supported systems in the presence of anisotropic orbits (see also Schwarzschild 1979, 1982; de Zeeuw 1985). In contrast, very little effort has been made in modeling rotation-dominated ellipticals, even though the entire low-mass end of the distribution of elliptical galaxies might be consistent with a picture of rotation-induced flattening (e.g., see Davies et al. 1983; and more recently Emsellem et al. 2011); similar comments apply to bulges.

For globular clusters, given the fact that they only exhibit modest amounts of flattening and given the success of the spherical King models, little work has been carried out in the direction of stationary self-consistent rotating models (with some notable exceptions, that is Wolley & Dickens 1962; Lynden-Bell 1962; Kormendy & Anand 1971; Lupton & Gunn 1987; and Lagoute & Longaretti 1996). Therefore, as far as rotation-dominated systems are concerned, much of the currently available modeling tools go back to the pioneering work of Prendergast & Tomer (1970), Wilson (1975), and Toomre (1982), intended to describe ellipticals, and of Jarvis & Freeman (1985) and Rowley (1988), devoted to bulges. In general, we may say that only very few rotating models with explicit distribution function are presently known (for a recent example, see Monari et al., in prep.). In this context, one should also mention the interesting work by Vandervoort (1980) on the collisionless analogues of the Maclaurin and Jacobi ellipsoids.

On the empirical side, a deeper study of quasi-relaxed rotating stellar systems is actually encouraged by the widespread conviction that the geometry of the inner parts of globular clusters, where some flattening is noted, should not be the result of tides, but rather of rotation (King 1961). In other words, it is frequently believed that tides and pressure anisotropy (and dust obscuration), even though playing some role in individual cases, should not be considered as the primary explanation of the observed flattening of Galactic globular clusters. Such conclusion is suggested by the White & Shawl (1987) database of ellipticities. However, for this database we note that: (i) the cluster flattening values do not all refer to a standard isophote, such as the cluster half-light radius (as also noted by van den Bergh 2008), (ii) the data mostly refer to the inner regions. Such limitations are crucial because there is observational evidence that the ellipticity of a cluster depends on radius (see Geyer et al. 1983). In turn, recent studies by Chen & Chen (2010), in which these restrictions are addressed more carefully, suggest that Galactic globular clusters are less round than previously reported, especially those in the region of the Galactic Bulge; at variance with previous analyses, it seems that their major axes preferentially point toward the Galactic Center, as naturally expected if their shape is of tidal origin.

The connection between flattening and internal rotation has been discussed in detail by means of nonspherical dynamical models in just a handful of cases, in particular for ω Cen (for an oblate rotator nonparametric model, see Merritt 1997; for an orbit-based analysis, see van de Ven 2006; for an application of the Wilson 1975, models; see Sollima et al. 2009), 47 Tuc (Meylan & Mayor 1986), M15 (van den Bosch et al. 2006), and M13 (Lupton et al. 1987). In addition, specifically designed 2D Fokker-Planck models (Fiestas et al. 2006) have been applied to the study of M5, NGC 2808, and NGC 5286. We recall that the detection of internal rotation in globular clusters is a challenging task, because the typical value of the ratio of the amplitude of the projected rotation velocity to the central velocity dispersion is only of a few tenths, for example V/σ₀ ≈ 0.46,0.32 for 47 Tuc and ω Cen, respectively (from a recent study by Bellazzini et al. 2012; for a summary of the results for several Galactic objects, see also Table 7.2 in Meylan & Heggie 1997). However, great progress made in the acquisition of photometric and kinematical information, and in particular of the proper motion of thousands of stars (for ω Cen, see van Leeuwen et al. 2000; Anderson & van der Marel 2010; for 47 Tucanae, see Anderson & King 2003; McLaughlin et al. 2006), makes this goal within reach (see Lane et al. 2009, 2010, for new kinematical measurements, in which rotation, when present, is clearly identified).

On the theoretical side, two general questions provide further motivation to study quasi-relaxed rotating stellar systems. On the one hand, many papers have studied the role of rotation in the general context of the dynamical evolution of globular clusters, but a solid interpretation is still missing. Early investigations (Agekian 1958; Shapiro & Marchant 1976) suggested that initially rotating systems should experience a loss of angular momentum induced by evaporation, that is, angular momentum would be removed by stars escaping from the cluster. Because of the small number of particles, N-body simulations were initially (Aarseth 1969; Wielen 1974; Akiyama & Sugimoto 1989) unable to clearly describe the complex interplay between relaxation and rotation. Later investigations, primarily based on a Fokker-Planck approach (Goodman 1983; Einsel & Spurzem 1999; Kim et al. 2002; Fiestas et al. 2006) have clarified this point, not only by testing the proposed mechanism of angular momentum removal by escaping stars, but also by showing that rotation accelerates the entire dynamical evolution of the system. More recent N-body simulations (Boily 2000; Ernst et al. 2007; Kim et al. 2008) confirm these conclusions and show that, when a three-dimensional tidal field is included, such acceleration is enhanced even further. The mechanism of angular momentum removal is generally considered to be the reason why Galactic globular clusters are much rounder than the (younger) clusters in the Magellanic Clouds, for which an age-ellipticity relation has been noted (Frenk & Fall 1982), but other mechanisms might operate to produce the observed correlations (Meylan & Heggie 1997; van den Bergh 2008).

On the other hand, the role of angular momentum during the initial stages of cluster formation should be better clarified. In the context of dissipationless collapse, relatively few investigations have considered the role of angular momentum in numerical experiments of violent relaxation (e.g., the pioneering studies by Gott 1973; see also Aguilar & Merritt 1990). Interestingly, the final equilibrium configurations resulting from such collisionless collapse show a central region with solid body rotation, while the external parts are characterized by differential rotation.

In conclusion, mastering the internal structure of spheroidal and triaxial stellar systems through a full spectrum of models, including rotation, is a prerequisite for studies of many empirical and theoretical issues. In addition, it is required for the interpretation of the relevant scaling laws (such as the Fundamental Plane, which appears to extend from the brightest, pressure supported ellipticals down to the low-luminosity end of the distribution of early-type galaxies, and possibly further down to the domain of globular clusters) and for investigations aimed at identifying the presence of invisible matter (in the form of central massive black holes or diffuse dark matter halos) from stellar dynamical measurements. It would thus be desirable to construct rotating models to be tested on low luminosity ellipticals, bulges, and globular clusters, especially now that important progress has been made in the collection and analysis of kinematical data. Presumably, many of these stellar systems are quasi-relaxed. In which directions should we explore deviations from the strictly relaxed case?

In the present paper, we consider two families of axisymmetric rotating models: the first one is characterized by the presence of solid-body rotation and isotropy in velocity space (see Appendix B in Paper I for a brief introduction). Indeed, full relaxation in the presence of nonvanishing total angular momentum suggests the establishment of solid-body rotation through the dependence of the relevant distribution function f = f(H) on the Jacobi integral H = E − ωJ_z (see Landau & Lifchitz 1967, p. 125). But, for applications to real stellar systems, one may take advantage of the fact that the collisional relaxation time may be large in the outer regions, so that in the outer parts the constraint of solid-body rotation might be released. In particular, for globular clusters we may argue that the outer parts fall into a tide-dominated regime, for which evaporation tends to erase systematic rotation even if initially present, as confirmed by the above-mentioned study based on the Fokker-Planck method. For the truncation, we may then consider a heuristic prescription to simulate the effects of tides, much like for the spherical King models.

In view of possible applications to globular clusters, we thus consider a second class of axisymmetric rotating models based on a distribution function dependent only on the energy and on the z-component of the angular momentum f = f(I) where I = I(E,J_z), with the property that I ~ E for stars with relatively high z-component of the angular momentum, while I ~ H = E − ωJ_z for relatively low values of J_z. Such models are indeed defined in order to have differential rotation, designed to be rigid in the center and to vanish in the outer parts, where the energy truncation becomes effective. As far as the velocity dispersion is concerned, this family may show a variety of profiles (depending on the values of the relevant free parameters), all of them characterized by the presence of isotropy in the central region. We thus add two classes of self-consistent models to the relatively short list of rotating stellar dynamical models currently available.

One aspect that plays an important role in defining a physically motivated distribution function, which often goes unnoticed (but see Hunter 1977; Davoust 1977; Rowley 1988), is the choice of the truncation prescription in phase space. The advantages and the limitations of alternative options available for the second family of models will be discussed in detail. In this context, we will also address the issue of whether these differentially rotating models fall within the class of systems for which rotation is constant on cylinders.

The paper is structured as follows. The properties of the family of rigidly rotating models, constructed on the basis of general statistical mechanics considerations, are illustrated in Sect. 2. The family of differentially rotating models, designed for the application to rotating globular clusters, is introduced in Sect. 3, where we briefly describe the method used for the solution of the self-consistent problem and discuss the relevant parameter space. Section 4 is devoted to a study of the intrinsic properties and Sect. 5 to the projected observables derived from differentially rotating models.

After illustrating in Sect. 6 the effect of different truncation prescriptions in phase space, we summarize the results and present our conclusions in Sect. 7. The appendices are devoted, respectively, to a discussion of the nonrotating limit of our families of models, to the details of the alternative truncation option for the second family, and to a description of the code used for the construction of the differentially rotating configurations. A study of the dynamical stability and of the long-term evolution of the families of models introduced here will be addressed by means of an extensive survey of specifically designed N-body simulations and will be presented in following papers (Varri et al. 2011; Varri et al., in prep.).

2. Rigidly rotating models

2.1. The distribution function

The construction of rigidly rotating configurations characterized by nonuniform density is a classical problem in the theory of rotating stars, starting with Milne (1923) and Chandrasekhar (1933), but it basically remained limited to the study of a fluid with polytropic equation of state, for which the solution of the relevant Poisson equation can be obtained by means of a semi-analytical approach (for a comprehensive description, see Chaps. 5 and 10 in Tassoul 1978; for an enlightening presentation of the general problem of rotating compressible masses, see Chap. 9 in Jeans 1929). The reader is referred to Paper I for a discussion of the application of some of the mathematical methods developed in that context to the construction of nonspherical truncated self-consistent stellar dynamical models. The deviations from spherical symmetry studied in Paper I are induced by the presence of a stationary perturbation characterized by a quadrupolar structure, that is, either an external tidal field or internal solid-body rotation. In particular, in Appendix B of Paper I we briefly outlined the extension of the King (1966) models to the case of internal solid-body rotation, of which, after a short summary of the relevant definitions, we now provide a full description in terms of the relevant intrinsic and projected properties.

The extension of the family of King models is performed by considering the distribution function: ${\mathit{f}}_{\mathrm{K}}^{\mathrm{r}}\mathrm{\left(}\mathit{H}\mathrm{\right)}\mathrm{=}\mathit{A}{\mathrm{e}}^{\mathrm{-}\mathit{a}{\mathit{H}}_{\mathrm{0}}}\mathrm{[}{\mathrm{e}}^{\mathrm{-}\mathit{a}\mathrm{(}\mathit{H}\mathrm{-}{\mathit{H}}_{\mathrm{0}}\mathrm{)}}\mathrm{-}{\mathrm{1}}^{\mathrm{]}}$(1)for H ≤ H₀ and ${\mathit{f}}_{\mathrm{K}}^{\mathrm{r}}\mathrm{\left(}\mathit{H}\mathrm{\right)}\mathrm{=}\mathrm{0}$ otherwise, where $\mathit{H}\mathrm{=}\mathit{E}\mathrm{-}\mathit{\omega }{\mathit{J}}_{\mathit{z}}$(2)denotes the Jacobi integral, with ω the angular velocity of the rigid rotation (hence the superscript r), assumed to take place around the z-axis. The quantities E and J_z are the specific one-star energy and z-component of the angular momentum, H₀ represents a cut-off constant of the Jacobi integral, while A and a are positive constants. In the corotating frame of reference the Jacobi integral can be written as the sum of the kinetic energy, the centrifugal potential Φ_cen(x,y) =  −(x² + y²)ω²/2 (with the equatorial plane (x,y) perpendicular to the rotation axis given by the z-axis), and the cluster mean-field potential Φ_C, to be determined self-consistently. Therefore, the dimensionless escape energy can be expressed as: $\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathit{a}\mathrm{\{}{\mathit{H}}_{\mathrm{0}}\mathrm{-}\mathrm{[}{\mathrm{\Phi }}_{\mathrm{C}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{cen}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{\left]}\mathrm{\right\}}\mathit{,}$(3)and the boundary of the cluster, implicitly defined by ψ(r) = 0, is an equipotential surface for the total potential Φ_C + Φ_cen. Note that, by construction, in the limit of vanishing internal rotation, this family of models reduces to the family of spherical King (1966) models (see Appendix A for a summary of the main properties of the family in the nonrotating limit).

The construction of the models requires the integration of the associated nonlinear Poisson equation, which, after scaling the spatial coordinates with respect to the scale length ${\mathit{r}}_{\mathrm{0}}\mathrm{=}{\left(\frac{\mathrm{9}}{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{0}}\mathit{a}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(4)can be written as $\mathrm{\nabla ̂}\mathrm{2}{\mathit{\psi }}^{\mathrm{\left(}\mathrm{int}\mathrm{\right)}}\mathrm{=}\mathrm{-}\mathrm{9}\left[\frac{\mathit{\rho ̂}\mathit{K}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{\left(}\mathrm{int}\mathrm{\right)}}\mathrm{\right)}}{\mathit{\rho ̂}\mathit{K}\mathrm{\left(}\mathrm{\Psi }\mathrm{\right)}}\mathrm{-}\mathrm{2}\mathit{\chi }\right]\mathit{,}$(5)where $\mathit{\chi }\mathrm{\equiv }\frac{{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{0}}}$(6)is the dimensionless parameter that characterizes the rotation strength and $\mathrm{\Psi }\mathrm{\equiv }{\mathit{\psi }}^{\mathrm{\left(}\mathrm{int}\mathrm{\right)}}\mathrm{\left(}\mathbf{0}\mathrm{\right)}$(7)is the depth of the potential well at the center. The dimensionless density profile is given by $\mathit{\rho ̂}\mathit{K}\mathrm{\left(}\mathit{\psi }\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathit{K}}\mathrm{\left(}\mathit{\psi }\mathrm{\right)}}{\mathit{Â}}\mathrm{=}{\mathrm{e}}^{\mathit{\psi }}\mathit{\gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}\mathit{,\psi }\right)\mathit{,}$(8)where $\mathit{Â}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }{\mathrm{2}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{A}}{\mathrm{3}{\mathit{a}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\mathit{a}{\mathit{H}}_{\mathrm{0}}}$(9)and γ denotes the incomplete gamma function. Therefore, the central density is given by . Outside the cluster, for negative values of ψ, we should refer to the Laplace equation $\mathrm{\nabla ̂}\mathrm{2}{\mathit{\psi }}^{\mathrm{\left(}\mathrm{ext}\mathrm{\right)}}\mathrm{=}\mathrm{18}\mathit{\chi }\mathit{.}$(10)The relevant boundary conditions are given by the requirement of regularity of the solution at the origin and by the condition that ψ^(ext) + aΦ_cen → aH₀ at large radii. The Poisson (internal) and Laplace (external) domains are thus separated by the boundary surface, which is unknown a priori. Therefore, we have to solve an elliptic partial differential equation in a free boundary problem. In particular, here we illustrate the properties of the solutions of the Poisson-Laplace equation obtained by using a perturbation method, which also requires an expansion of the solution in Legendre series. To obtain a uniformly valid solution over the entire space, an asymptotic matching is performed between the internal and the external solution, using the Van Dyke principle (see Van Dyke 1975). This method of solution is basically the same as proposed by Smith (1975) for the construction of rotating configurations with polytropic index n = 3/2. We have calculated the complete solution up to second-order in the rotation strength parameter χ. The final solution is expressed in spherical coordinates and the resulting configurations are characterized by axisymmetry (i.e., the density distribution and the potential do not depend on the azimuthal angle φ). For the details of the method for the construction of the solution the reader is referred to Paper I and to Varri (2012), in which the complete calculation is provided.

2.2. The parameter space

The resulting models are characterized by two dimensional scales (e.g., the total mass and the core radius) and two dimensionless parameters. As in the spherical King models, the first parameter measures the concentration of the configuration; we thus consider the quantity¹ Ψ (see Eq. (7)) or equivalently c ≡ log (r_tr/r₀), where r_tr is the truncation radius of the spherical King model associated with a given value of Ψ. The second dimensionless parameter χ (see Eq. (6)) characterizes the rotation strength measured in terms of the frequency associated with the central density of the cluster.

For every value of the dimensionless central concentration Ψ there exists a maximum value of the rotation strength parameter, corresponding to a critical model, for which the boundary is given by the critical constant-ψ surface. The boundary surface of a representative critical model (with Ψ = 2) is depicted in Fig. 1; the surface is such that all the points on the equatorial plane are saddle points, where the centrifugal force balances the self-gravity. We refer to their distance from the origin as , the break-off radius. In the constant-ψ family of surfaces associated with a given value of Ψ, the critical surface thus separates the open from the closed surfaces. Consistent with the assumption of stationarity, only configurations bounded by closed surfaces are considered here. This unique geometrical characterization suggests that the effect of the rotation may be expressed also in terms of an extension parameter $\mathit{\delta }\mathrm{=}\frac{\mathit{r̂}\mathrm{tr}}{\mathit{r̂}\mathrm{B}}\mathit{,}$(11)which provides an indirect measure of the deviations from sphericity of a configuration, by considering the ratio between the truncation radius of the corresponding spherical King model and the break-off radius of the associated critical surface. Therefore, a given model may be labelled by the pair (Ψ,χ) or equivalently by the pair (Ψ,δ). For a given Ψ, there is thus a maximum value of the allowed rotation, which we may express as χ_cr or δ_cr. A model with δ < δ_cr may be called subcritical.

For each value of Ψ, the critical value of the rotation parameter can be found by numerically solving the system²$\{\begin{array}{c}\\ {\mathit{\partial }}_{\mathit{r̂}}\hspace{0.17em}\mathit{\psi }\mathrm{(}\mathit{r̂}\mathrm{=}\mathit{r̂}\mathrm{B}\mathit{,\theta }\mathrm{=}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{;}{\mathit{\chi }}_{\mathrm{cr}}\mathrm{)}\mathrm{=}\mathrm{0}\\ \mathit{\psi }\mathrm{(}\mathit{r̂}\mathrm{=}\mathit{r̂}\mathrm{B}\mathit{,\theta }\mathrm{=}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{;}{\mathit{\chi }}_{\mathrm{cr}}\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(12)where the unknowns are and χ_cr. In terms of the extension parameter, for a given Ψ, the critical condition occurs when . This value is obtained by inserting in Eq. (12) the zeroth-order expression for the cluster potential, as discussed in detail in Sect. 2 of Paper II for the tidal problem (see pp. 251–255 in Jeans 1929, for an equivalent discussion referred to the purely rotating Roche model, that is a rotating configuration in which a small region with infinite density is surrounded by an “atmosphere” of negligible mass).

The parameter space for the second-order models is presented in Fig. 2 (which corresponds to Fig. 1 of Paper II describing the tidal models). Two rotation regimes exist, namely the regime of low-deformation (δ ≪ δ_cr, bottom left corner), where internal rotation does not affect significantly the morphology of the configuration, which remains very close to spherical symmetry, and that of high-deformation (δ ≈ δ_cr, close to the solid line), where the model is highly affected by the nearly critical rotation velocity, especially in the outer parts. Note that the actual regime depends on the combined effect of rotation strength and of concentration. In other words, the models described here belong to the class of rotating configurations characterized by equatorial break-off (“region of equatorial break-off”, see Fig. 44, p. 267 in Jeans 1929), for which the limiting case is given by the purely rotating Roche model.

A global kinematical characterization, complementary to the information provided by the rotation strength parameter χ, is offered by the parameter t = K_ord/ | W | , defined as the ratio between ordered kinetic energy and gravitational energy. Figure 3 illustrates the relation between the two parameters for models with selected values of Ψ and increasing values of χ, up to the critical configuration characterized by χ_cr. The parameter t increases linearly for increasing values of χ (with a slope dependent on the concentration parameter Ψ). A similar linear relation is observed for small values of eccentricity (e ≪ 1) in the sequence of Maclaurin oblate spheroids t(e) ~ χ(e) ~ 2e²/15 (for the definitions of the two parameters in the context of Maclaurin spheroids, see Eqs. (10.20) and (10.24) in Bertin 2000); in Fig. 3 the relevant linear relation is normalized with respect to the maximum value of the rotation strength parameter attained in the sequence of Maclaurin spheroids χ_max = 0.11233 (see Eq. (10), p. 80 in Chandrasekhar 1969).

2.3. Intrinsic properties

The geometry of the models, reflecting the properties of the centrifugal potential, is characterized by symmetry around the ẑ-axis and reflection symmetry with respect to the equatorial plane . As expected, compared to the corresponding spherical King models, the rotating models stretch out on the equatorial plane and are slightly flattened along the direction of the rotation axis (see Fig. 4 for the density profiles of selected critical second order models evaluated on the equatorial plane and along the ẑ-axis). In general, configurations in the low-deformation regime (δ ≪ δ_cr), regardless of the value of concentration Ψ, are almost indistinguishable from the corresponding spherical King models.

Models in the intermediate and high-deformation regime (δ ≈ δ_cr) show modest deviations from spherical symmetry in the central regions, while they are significantly flattened in the outer parts. The intrinsic eccentricity profile, defined as , where â and are semi-major and semi-minor axes of the isodensity surfaces, is a monotonically increasing function of the semi-major axis (see Fig. 5 for the eccentricity profiles of selected critical second-order models). We recall that the geometry of the boundary surface of a critical model depends only slightly on the value of the concentration parameter. In particular, in the critical case, the break-off radius represents the distance from the center of the outermost points of the boundary surface on the equatorial plane and the truncation radius is approximately the distance of the last point on the polar axis (i.e., the ẑ-axis). Therefore, since , the value of the termination points of the eccentricity profiles of critical models is approximately the same (e ≈ 0.75; see the termination points of the solid lines in Fig. 5).

In addition, by using the multipolar structure of the solution of the Poisson-Laplace equation obtained with the perturbation method, the asymptotic behavior of the eccentricity profiles in the central regions can also be evaluated analytically. Since the distribution function depends only on the Jacobi integral (i.e., the isolating energy integral in the corotating frame of reference), the density and the velocity dispersion profiles are functions of only the escape energy (see Eqs. (8) and (17), respectively). Therefore, there is a one-to-one correspondence between equipotential, isodensity, and isobaric surfaces and the eccentricity profiles can be calculated with reference to just one of these families of surfaces. In fact, if expanded to second order in the dimensionless radius, the escape energy in the internal region reduces to: $\begin{array}{ccc}{\mathit{\psi }}^{\mathrm{\left(}\mathrm{int}\mathrm{\right)}}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}& \mathrm{=}& \mathrm{\Psi }\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}[\mathrm{-}\mathrm{3}\mathrm{+}\mathrm{6}\mathit{\chi }\mathrm{2}{\mathit{A}}_{\mathrm{2}}{\mathit{U}}_{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{\chi }\\ & & \mathrm{+}\hspace{0.17em}\mathrm{(}{\mathit{B}}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{U}}_{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{\chi }}^{\mathrm{2}}\mathrm{]}\mathit{r̂}\mathrm{2}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{4}\mathrm{\right)}\mathit{,}\end{array}$(13)where U₂(θ) denotes the normalized Legendre polynomial with l = 2 and A₂,B₂ are appropriate (negative) coefficients, depending on Ψ, which are determined by asymptotically matching the internal and external solution, in order to have continuity on the entire domain (see Eqs. (62) and (68) in Paper I, to be interpreted as indicated in Appendix B of Paper I). By setting , we thus find that, in the innermost region, the eccentricity tends to the following nonvanishing central value ${\mathit{e}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{[}\mathrm{6}{\mathit{A}}_{\mathrm{2}}\hspace{0.17em}\mathit{\chi }\mathrm{+}\mathrm{3}\mathrm{(}{\mathit{B}}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{\chi }}^{\mathrm{2}}{\mathrm{]}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{\left[}\mathrm{6}\sqrt{\mathrm{2}\mathit{/}\mathrm{5}}\mathrm{\right(}\mathrm{2}\mathit{\chi }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{+}\mathrm{4}{\mathit{A}}_{\mathrm{2}}\hspace{0.17em}\mathit{\chi }\mathrm{+}\mathrm{2}\mathrm{(}{\mathit{B}}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{\chi }}^{\mathrm{2}}{\mathrm{]}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}$(14)Therefore, the central value of the eccentricity is finite, of order , and strictly vanishes only in the limit of vanishing rotation strength. This result is nontrivial because the centrifugal potential (which induces the deviations from sphericity) is a homogeneous function of the spatial coordinates. Therefore, we might naively expect that, in their central regions, the models reduce to a perfectly spherical shape (i.e., e₀ = 0), even for finite values of the rotation strength. This property has been noted also in the family of triaxial tidal models, in which the tidal potential plays the role of the centrifugal potential (see Sect. 3.1 of Paper II).

Deviations from spherical symmetry can also be described, in a global way, by the quadrupole moment tensor, defined as $\begin{array}{ccc}{\mathit{Q}}_{\mathit{ij}}& \mathrm{=}& \int \mathit{V}\hspace{0.17em}\mathrm{(}\mathrm{3}{\mathit{x}}_{\mathit{i}}{\mathit{x}}_{\mathit{j}}\mathrm{-}{\mathit{r}}^{\mathrm{2}}{\mathit{\delta }}_{\mathit{ij}}\mathrm{)}\mathit{\rho }\mathrm{\left(}{r}\mathrm{\right)}{\mathrm{d}}^{\mathrm{3}}\mathit{r}\\ & \mathrm{=}& \mathit{Â}{\mathit{r}}_{\mathrm{0}}^{\mathrm{5}}{\mathrm{\int }}_{\mathit{V}}\hspace{0.17em}\mathrm{(}\mathrm{3}\mathit{x̂}\mathit{i}\mathit{x̂}\mathit{j}\mathrm{-}\mathit{r̂}\mathrm{2}{\mathit{\delta }}_{\mathit{ij}}\mathrm{)}\mathit{\rho ̂}\mathrm{\left(}{r̂}\mathrm{\right)}{\mathrm{d}}^{\mathrm{3}}\mathit{r̂}\mathrm{=}\mathit{Â}{\mathit{r}}_{\mathrm{0}}^{\mathrm{5}}\mathit{Q̂}\mathit{ij}\mathit{,}\end{array}$(15)where the integration is performed in the volume V of the entire configuration. It can be easily shown that in our coordinate system the tensor is diagonal and that Q_xx = Q_yy. The components of the tensor can be calculated explicitly $\mathit{Q̂}\begin{array}{}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ \mathit{xx}\end{array}\mathrm{=}\mathit{Q̂}\begin{array}{}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ \mathit{yy}\end{array}\mathrm{=}\mathrm{-}\frac{\mathit{Q̂}\begin{array}{}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ \mathit{zz}\end{array}}{\mathrm{2}}\mathrm{=}\mathrm{2}\mathit{\pi }\frac{\sqrt{\mathrm{1}}\mathrm{0}}{\mathrm{9}}\mathit{\rho ̂}\mathit{K}\mathrm{\left(}\mathrm{\Psi }\mathrm{\right)}\left({\mathit{a}}_{\mathrm{2}}\mathit{\chi }\mathrm{+}{\mathit{b}}_{\mathrm{2}}\frac{{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right)\mathit{,}$(16)the quantities a₂ and b₂ are appropriate (positive) coefficients, depending on Ψ, resulting from the asymptotic matching of the internal and external solution of the Poisson-Laplace equation. The sign of the components are consistent with the above mentioned compression and stretching of the density distribution. The expression in Eq. (16) is calculated from the second-order external solution; the first-order expression is recovered by dropping the quadratic term in the parameter χ. At variance with the tidal case, the ratio is independent of the rotation parameter χ. This analytical result has been compared to the ratio of the components of the quadrupole tensor determined by direct numerical integration (performed by means of the algorithm VEGAS, see Press et al. 1992) and good agreement has been found³. For the tidal case, the detailed calculation can be found in Sect. 3.3 and Appendix B of Paper II; such calculation is easily adapted to the rotating case.

By construction, the models are isotropic in velocity space, with the dimensionless scalar velocity dispersion given by $\mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathit{K}\end{array}\mathrm{\left(}\mathit{\psi }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{\mathrm{5}}\frac{\mathit{\gamma }\mathrm{(}\mathrm{7}\mathit{/}\mathrm{2}\mathit{,\psi }\mathrm{)}}{\mathit{\gamma }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}\mathit{,\psi }\mathrm{)}}\mathrm{·}$(17)As noted for the intrinsic density profiles, a compression along the vertical axis and a stretching along the equatorial plane occur also for the velocity dispersion profiles (the profiles of selected critical second-order models are shown in Fig. 6).

The mean rotation velocity (which is subtracted away when the corotating frame of reference is considered) characterizing the models is defined as $\mathrm{⟨}{v}\mathrm{⟩}\mathrm{=}\mathit{\omega }\hspace{0.17em}\hat{{e}}{}_{\mathit{z}}\mathrm{\times }\mathit{r}\hspace{0.17em}\hat{{e}}{}_{\mathit{r}}\mathrm{=}\mathrm{⟨}{\mathit{v}}_{\mathit{\phi }}\mathrm{⟩}\hat{{e}}{}_{\mathit{\phi }}$; in the adopted dimensionless units, the azimuthal component can be written as $\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}\mathrm{=}\mathrm{⟨}{\mathit{v}}_{\mathit{\phi }}\mathrm{⟩}{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{3}{\mathit{\chi }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{r̂}\sin \mathit{\theta }\mathit{.}$(18)As expected, the mean velocity is constant on cylinders (in our coordinate system, the cylindrical radius is defined by ). The relevant dimensionless angular velocity is linked to the rotation strength parameter by the following relation (the numerical factor 3 is due to the adopted scale length, see Eq. (4)). The rotation profiles of selected second-order critical models are represented in Fig. 7. For all the models, as we approach the boundary of the configuration, the ratio quickly diverges since at the boundary the rotation velocity tends to a finite value while the velocity dispersion vanishes; this behavior is observed in every direction (except for the ẑ-axis, on which the rotation velocity vanishes by definition).

2.4. Projected properties

For a comparison of the models with the observations (under the assumption of a constant mass-to-light ratio), we have then computed surface (projected) density profiles and isophotes. The projection has been performed along selected directions, identified by the viewing angle (θ,φ) corresponding to the ẑ_P axis of a new coordinate system related to the intrinsic system by the transformation ; the rotation matrix R = R₁(θ)R₃(φ) is expressed in terms of the viewing angles, by taking the axis as the line of nodes (i.e., the same projection rule we adopted for the triaxial tidal models, see Sect. 3.4 in Paper II). Since the rigidly rotating models are characterized by axisymmetry with respect to the ẑ-axis and reflection symmetry with respect to the equatorial plane, it is sufficient to choose the viewing angles from the -plane of the intrinsic coordinate system. In particular, we used the line of sights defined by θ_i = i(π/8) and φ = 0 with i = 0,...,4, and we calculated (by numerical integration, using the Romberg’s rule) the dimensionless projected density $\mathrm{\Sigma ̂}\mathrm{\left(}\mathit{x̂}\mathrm{P}\mathit{,}\mathit{ŷ}\mathrm{P}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathit{ẑ}\mathrm{sp}}^{\mathit{ẑ}\mathrm{sp}}\mathit{\rho ̂}\mathrm{\left(}{r̂}\mathrm{P}\mathrm{\right)}\mathrm{d}\mathit{ẑ}\mathrm{P}\mathit{,}$(19)where $\mathit{ẑ}{}_{\mathrm{sp}}\mathrm{=}\mathrm{(}\mathit{x̂}\begin{array}{}\mathrm{2}\\ \mathrm{e}\end{array}\mathrm{-}\mathit{x̂}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}\mathrm{-}\mathit{ŷ}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ with the edge of the model along the axis of the intrinsic coordinate system. The projection plane has been sampled on an equally-spaced square cartesian grid centered at the origin.

The first five panels of Fig. 8 show the projected images of a critical second-order model with Ψ = 2; the first and the fifth panel correspond, respectively, to the least and to the most favorable line of sight for the detection of the intrinsic flattening of the model (the ẑ and -axis of the intrinsic coordinate system, that is “face-on” and “edge-on” view).

The morphology of the isophotes of a given projected image can be described in terms of the ellipticity profile, defined as where â_P and are the principal semi-axes, as a function of the semi-major axis â_P. As already noted for the (intrinsic) eccentricity profile, the deviation from circularity increases with the distance from the origin. In the inner region, the central value of the ellipticity is consistent with the central eccentricity e₀ calculated in the previous subsection. The last panel of Fig. 8 illustrates the ellipticity profiles corresponding to the projections displayed in the previous panels. In addition, the isophotes of models in the high deformation regime (δ ≈ δ_cr), if projected along appropriate line of sights, show clear departures from a pure ellipse, that can be characterized as a “disky” overall trend (e.g., see Jedrzejewski 1987), particularly evident in the outer parts (see fourth and fifth panels in Fig. 8).

Because the family of rigidly rotating models is characterized by simple kinematical properties (pressure isotropy and solid-body rotation), that have been already presented in detail with reference to three-dimensional configurations, for brevity, the derivation of the projected kinematical properties is omitted here; a full description can be found in Varri (2012).

3. Differentially rotating models

3.1. Choice of the distribution function

As indicated in the Introduction, theoretical and observational motivations have brought us to look for more realistic configurations, characterized by differential rotation. Thus we focus our attention on axisymmetric systems, within the class of distribution functions that depend only on the energy E and the z-component of the angular momentum J_z, and we consider the integral $\mathit{I}\mathrm{\left(}\mathit{E,}{\mathit{J}}_{\mathit{z}}\mathrm{\right)}\mathrm{=}\mathit{E}\mathrm{-}\frac{\mathit{\omega }{\mathit{J}}_{\mathit{z}}}{\mathrm{1}\mathrm{+}\mathit{b}{\mathit{J}}_{\mathit{z}}^{\mathrm{2}\hspace{0.17em}\mathit{c}}}\mathit{,}$(20)where ω, b, and c > 1/2 are positive constants. The quantity I(E,J_z) reduces to the Jacobi integral for small values of the z-component of the angular momentum and tends to the single-star energy in the limit of high values of J_z. Therefore, if we refer to a distribution function of the form f = f(I), we may argue that ω is related to the angular velocity in the central region of the system, characterized by approximately solid-body rotation, whereas the positive constants b,c will determine the shape of the radial profile of the rotation profile. In view of the arguments that have led to the truncation prescription that characterizes King model, we decided to introduce a truncation in phase-space based exclusively on the single-star energy with respect to a cut-off constant E₀.

For simplicity, we consider two families of distribution functions. The first family is defined as ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}\mathrm{=}\mathit{A}{\mathrm{e}}^{\mathrm{-}\mathit{a}{\mathit{E}}_{\mathrm{0}}}\mathrm{[}{\mathrm{e}}^{\mathrm{-}\mathit{a}\mathrm{(}\mathit{I}\mathrm{-}{\mathit{E}}_{\mathrm{0}}\mathrm{)}}\mathrm{-}\mathrm{1}\mathrm{+}\mathit{a}\mathrm{(}\mathit{I}\mathrm{-}{\mathit{E}}_{\mathrm{0}}{\mathrm{\left)}}^{\mathrm{\right]}}$(21)if E ≤ E₀ and ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}\mathrm{=}\mathrm{0}$ otherwise, so that both ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ and its derivative with respect to E are continuous. We refer to this truncation prescription as Wilson truncation (hence the subscript WT) because, in the limit of vanishing internal rotation (ω → 0), this family reduces to the spherical limit of the distribution function proposed by Wilson (1975). The superscript d in Eq. (21) indicates the presence of differential rotation.

The second family is defined by the distribution function ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}\mathrm{=}\mathit{A}{\mathrm{e}}^{\mathrm{-}\mathit{a}{\mathit{E}}_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{-}\mathit{a}\mathrm{(}\mathit{I}\mathrm{-}{\mathit{E}}_{\mathrm{0}}\mathrm{)}}$(22)if E ≤ E₀ and ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}\mathrm{=}\mathrm{0}$ otherwise; therefore the function, characterized by plain truncation (hence the subscript PT), is discontinuous with respect to E. In the limit of vanishing internal rotation, it reduces to the spherical limit of the function proposed by Prendergast & Tomer (1970), which leads to the truncated isothermal sphere (see also Woolley & Dickens 1962). A summary of the main properties and definitions of the relevant nonrotating limit of the two families of models is presented in Appendix A.

In both cases the distribution functions are positive definite ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}\mathit{,}\hspace{0.17em}{\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}\mathrm{\ge }\mathrm{0}$, by construction. Curiously, a naive extension of King models $\mathit{f}\mathrm{=}{\mathit{f}}_{\mathrm{K}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ with a similar truncation in energy alone would define a distribution function that is not positive definite in the whole domain of definition.

Sharp gradients or discontinuities in phase space (such as the ones associated with the truncation prescription of ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$) are expected to be associated with evolutionary processes dictated either by collective modes or by any small amount of collisionality. Therefore, the first truncation prescription, corresponding to a smoother distribution in phase space, is to be preferred from a physical point of view as the basis for a realistic equilibrium configuration (in principle, we might have referred to even smoother functions; see Davoust 1977). In addition, a full analysis of the configurations defined by ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ shows that this family of models exhibits a number of interesting intrinsic and projected properties, more appropriate for application to globular clusters, with respect to the models defined by ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$.

Therefore, the following Sects. 4 and 5 are devoted to the full characterization of the family of models defined by ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ (for a summary of the properties of the families of models studied in the present paper, see Table 1). The intrinsic properties of the family of models defined by ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ are summarized in Appendix B. In this investigation, we decided to briefly mention and to keep also the second family not only because it extends a well-known family of models, but also because it allows us to check directly an important aspect of model construction that had been noted by Hunter (1977). This is that the truncation prescription affects the density distribution in the outer parts of the models significantly. This general point is even more important if we are interested in modeling the outermost regions of globular clusters (with particular reference to the studies of the “extra-tidal lights”, e.g., see Jordi & Grebel 2010).

3.2. The construction of the models

The construction of the models requires the integration of the relevant Poisson equation, supplemented by a set of boundary conditions equivalent to the one described in Sect. 2 for rigidly rotating models. In this case, we obtain the solution by means of an iterative approach, based on the method proposed by Prendergast & Tomer (1970), in which an improved solution ψ⁽ⁿ⁾ of the dimensionless Poisson equation is obtained by evaluating the source term on the right-hand side with the solution from the immediately previous step (for an application of the same method to the construction of configurations shaped by an external tidal field, see Sect. 5.2 in Paper I) $\mathrm{\nabla ̂}\mathrm{2}{\mathit{\psi }}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{=}\mathrm{-}\frac{\mathrm{9}}{\mathit{\rho ̂}\mathrm{0}}\mathit{\rho ̂}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,}{\mathit{\psi }}^{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{\right)}\mathrm{;}$(23)here the dimensionless escape energy is given by $\mathit{\psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathit{a}\mathrm{[}{\mathit{E}}_{\mathrm{0}}\mathrm{-}{\mathrm{\Phi }}_{\mathit{C}}\mathrm{(}{r}\mathrm{\left)}\mathrm{\right]}\mathit{,}$(24)the dimensionless radius is defined as , with the same scale length introduced in Eq. (4), and indicates the dimensionless central density. The relevant density profile , with  defined as in Eq. (9), results from the integration in velocity space of the distributions function defined by ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$. It is clear that the general strategy for the construction of the self-consistent solution is applicable also to the density derived from ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$.

The iteration is seeded by the corresponding spherical models, that is, the Wilson and the Prendergast-Tomer spheres respectively, and is stopped when numerical convergence is reached (see Appendix C for details). At each iteration step, the scheme requires the expansion in Legendre series of the density and the potential $\begin{array}{ccc}{\mathit{\psi }}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}{r̂}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathit{\psi }}_{\mathit{l}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}{\mathit{U}}_{\mathit{l}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}\mathit{,}\\ \mathit{\rho ̂}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{\left(}{r̂}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\mathit{\rho ̂}\begin{array}{}\mathrm{\left(}\mathit{n}\mathrm{\right)}\\ \mathit{l}\end{array}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}{\mathit{U}}_{\mathit{l}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}\mathit{.}\end{array}$The associated Cauchy problems for the radial functions ${\mathit{\psi }}_{\mathit{l}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}$ are therefore $\left[\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}\mathit{r̂}\mathrm{2}}\mathrm{+}\frac{\mathrm{2}}{\mathit{r̂}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r̂}}\mathrm{-}\frac{\mathit{l}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathit{r̂}\mathrm{2}}\right]{\mathit{\psi }}_{\mathit{l}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{=}\mathrm{-}\frac{\mathrm{9}}{\mathit{\rho ̂}\mathrm{0}}\mathit{\rho ̂}\begin{array}{}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\\ \mathit{l}\end{array}\mathit{,}$(27)supplemented by the following boundary conditions $\begin{array}{ccc}{\mathit{\psi }}_{\mathrm{0}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}& \mathrm{=}& \mathrm{\Psi }\sqrt{\mathrm{2}}\mathit{,}\\ {\mathit{\psi }}_{\mathit{l}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}& \mathrm{=}& \mathrm{0}\mathit{,}\mathrm{for}\mathit{l}\ne \mathrm{0}\\ {{\mathit{\psi }}_{\mathrm{0}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}& \mathrm{=}& {{\mathit{\psi }}_{\mathit{l}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$

where Ψ is the depth of the dimensionless potential well at the center. By using the method of variation of arbitrary constants, the radial functions can be expressed in integral form as follows $\begin{array}{ccc}{\mathit{\psi }}_{\mathrm{0}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}& \mathrm{=}& \mathrm{\Psi }\sqrt{\mathrm{2}}\mathrm{-}\frac{\mathrm{9}}{\mathit{\rho ̂}\mathrm{0}}[{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r̂}}\mathit{r̂}\mathrm{\prime }\mathit{\rho ̂}\begin{array}{}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\\ \mathrm{0}\end{array}\mathrm{\left(}\mathit{r̂}\mathrm{\prime }\mathrm{\right)}\mathrm{d}\mathit{r̂}\mathrm{\prime }\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathit{r̂}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r̂}}\mathit{r̂}\mathrm{\prime }\mathrm{2}\mathit{\rho ̂}\begin{array}{}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\\ \mathrm{0}\end{array}\mathrm{\left(}\mathit{r̂}\mathrm{\prime }\mathrm{\right)}\mathrm{d}\mathit{r̂}\mathrm{\prime }]\mathit{,}\\ {\mathit{\psi }}_{\mathit{l}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{9}}{\mathrm{(}\mathrm{2}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{\rho ̂}\mathrm{0}}[\mathit{r̂}\mathit{l}{\mathrm{\int }}_{\mathit{r̂}}^{\mathrm{\infty }}\mathit{r̂}\mathrm{\prime }\mathrm{1}\mathrm{-}\mathit{l}\mathit{\rho ̂}\begin{array}{}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\\ \mathit{l}\end{array}\mathrm{\left(}\mathit{r̂}\mathrm{\prime }\mathrm{\right)}\mathrm{d}\mathit{r̂}\mathrm{\prime }\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathit{r̂}\mathit{l}\mathrm{+}\mathrm{1}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r̂}}\mathit{r̂}\mathrm{\prime }\mathit{l}\mathrm{+}\mathrm{2}\mathit{\rho ̂}\begin{array}{}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\\ \mathit{l}\end{array}\mathrm{\left(}\mathit{r̂}\mathrm{\prime }\mathrm{\right)}\mathrm{d}\mathit{r̂}\mathrm{\prime }]\mathit{.}\end{array}$The factor $\sqrt{\mathrm{2}}$ appearing in Eqs. (28) and (31) is due to the normalization assumed for the Legendre polynomials.

3.3. The parameter space

Much like in the case of rigidly rotating models, in both families ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ and ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$, the resulting models are characterized by two scales, associated with the positive constants A and a, and two dimensionless parameters (Ψ,χ), measuring concentration and central rotation strength, respectively. In addition, two new dimensionless parameters, namely c and $\mathit{b̅}\mathrm{=}\mathit{b}\hspace{0.17em}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}\mathit{c}}\hspace{0.17em}{\mathit{a}}^{\mathrm{-}\mathit{c}}\mathit{,}$(33)determine the shape of the rotation profile. Variations in the parameter and c are found to be less important. Minor changes in the model properties are found up to , above which the precise value of c has only little impact on the properties of the configurations.

For each family of differentially rotating models, for given values there exists a maximum value of the central rotation strength parameter χ_max, corresponding to the last configuration for which the iteration described in Sect. 3.2 converges. Such maximally rotating configurations exhibit highly deformed morphologies, characterized by the presence of a sizable central toroidal structure, which will be described in detail in the following sections.

As for the parameter space of rigidly rotating models, it is useful to introduce different rotation regimes, defined on the basis of the deviations from spherical symmetry introduced by the presence of differential rotation. With particular reference to the parameter space of the models defined by ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$, we introduce some threshold values in central dimensionless angular velocity, which, in this family of models, is related to the central rotation strength parameter χ by the relation , as for the rigidly rotating models. In particular, configurations in the moderate rotation regime have (the thin-striped area in Fig. 9), are quasi-spherical in the outer parts, while they are progressively more flattened when approaching the central region, as the value of χ increases. For the models falling in this rotation regime, the central toroidal structure is absent or, when low values of the concentration parameter Ψ are considered, not significant. Configurations with (the wide-striped area in Fig. 9) are defined as rapidly rotating models. The extreme rotation regime is defined by the condition (the gray area in Fig. 9); in this case, the models always show a central toroidal structure, which becomes more extended as the central rotation strength increases. In particular, in the last regime, the entire volume of a configuration is dominated by the central toroidal structure.

As for the rigidly rotating models described in Sect. 2, a global kinematical characterization is offered by the parameter t = K_ord/ | W | , defined as the ratio between ordered kinetic energy and gravitational energy. Figure 10 illustrates the relation between the two parameters, for models with selected values of Ψ, , and c. Note that the transition from rapid to extreme rotation corresponds to values of the parameter t in the range  [0.075,0.135]  (the precise value depends on the value of Ψ). A naive application of the Ostriker & Peebles (1973) criterion, which states that axisymmetric stellar systems with t > 0.14 are dynamically unstable with respect to bar modes, would suggest that the majority of the models in the extreme rotation regime are dynamically unstable. A detailed stability analysis of the configurations in the three rotation regimes has been performed by means of specifically designed N-body simulations and will be presented in a separated paper (Varri et al. 2011, and in prep.).

To provide a systematic description of the intrinsic and projected properties, we will study the equilibrium configurations as sequences of models characterized by a given value of concentration Ψ, in the range  [1,7] , and increasing values of χ, up to the maximum value allowed; such sequences are constructed by fixing (unless otherwise stated).

4. Intrinsic properties of the differentially rotating models

4.1. The intrinsic density profile

The relevant density profile is obtained by integration in velocity space of the distribution function ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ (see Eq. (21)). It is convenient to introduce in the velocity space a spherical coordinate system (v,μ,λ), in which v is magnitude of the velocity vector, while μ and λ are the polar and azimuthal angle, respectively. After some manipulation, the density profile can be expressed in dimensionless form as $\begin{array}{ccc}\mathit{\rho ̂}\mathrm{WT}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{3}}{\mathrm{4}}{\mathrm{e}}^{\mathit{\psi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\psi }}\mathrm{d}\mathit{s}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{s}}{\mathit{s}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathrm{d}\mathit{tg}\mathrm{\left(}\mathit{s,t,}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\\ & & \mathrm{-}{\mathit{\psi }}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(34)where the function in the integrand is defined as $\mathit{g}\mathrm{\left(}\mathit{s,t,}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\exp \left(\frac{\mathrm{3}{\mathit{\chi }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{t}\hspace{0.17em}\mathit{r̂}\sin \mathit{\theta }\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{s}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{1}\mathrm{+}\mathit{b̅}{\left[\mathit{t}\hspace{0.17em}\mathit{r̂}\sin \mathit{\theta }\hspace{0.17em}\mathrm{(}\mathrm{2}\mathit{s}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\right]}^{\mathrm{2}\mathit{c}}}\right)\mathrm{;}$(35)for completeness, we note that the two dimensionless variables in the double integral can be expressed in terms of the previous variables as t = cosμ and s = av²/2. Because the distribution function depends only on the energy E and the z-component of the angular momentum J_z, the resulting models are axisymmetric and therefore the density profile depends only on the radius and the polar angle θ. The density depends on the spatial coordinates explicitly and implicitly, through the dimensionless escape energy ; such explicit dependence is the reason why, in this case, the isodensity, equipotential, and isobaric surfaces are not in one-to-one correspondence, at variance with the family of rigidly rotating models. The presence in Eq. (34) of the terms with fractional powers of ψ is due to the adopted truncation in phase space; in particular, it is directly related to the presence of the terms e ^− aE₀ [−1 + a(I − E₀)]  in Eq. (21).

The central value of the density profile depends only on the concentration parameter Ψ and is given by , consistent with the central value of the density profile obtained in the nonrotating limit (see Eq. (A.4) in Appendix A). This result corresponds to the fact that the integral I(E,J_z) (see Eq. (20)) reduces to the Jacobi integral for small values of J_z, which implies that the rotation is approximately rigid in the central regions of a configuration (see the next subsection for details); therefore, the mean rotation velocity vanishes at the origin and the density distribution reduces to its nonrotating limit.

The double integral in Eq. (34) requires a numerical integration (see Appendix C for details). Some insight into the behavior of the density profile can be gained by calculating the relevant asymptotic expansions in the central region and in the outer parts. Around the origin, up to second order in radius, the density profile reduces to $\begin{array}{ccc}\mathit{\rho ̂}\mathrm{WT}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,}\mathrm{\Psi }\mathrm{\right)}& \mathrm{=}& \mathit{\rho ̂}\mathrm{WT}\mathit{,}\hspace{0.17em}\mathrm{0}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{\Psi }}\mathit{\gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}\mathit{,}\mathrm{\Psi }\right)\mathrm{[}\mathrm{9}\mathit{\chi }{\sin }^{\mathrm{2}}{\mathit{\theta }}^{}\\ & & \mathrm{+}{\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }\mathit{r̂}\mathrm{2}}|}_{\mathrm{0}}]\mathit{r̂}\mathrm{2}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{4}\mathrm{\right)}\mathit{,}\end{array}$(36)which depends explicitly on the concentration Ψ and the central rotation strength χ, and implicitly on the parameters and c, through the second order derivative of the escape energy evaluated at . Such derivative can be calculated from the radial functions given in Eqs. (31), (32) $\begin{array}{ccc}{\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }\mathit{r̂}\mathrm{2}}|}_{\mathrm{0}}& \mathrm{=}& {\mathit{\psi }}_{\mathrm{0}}^{\mathrm{\prime \prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathit{U}}_{\mathrm{0}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{+}{\mathit{\psi }}_{\mathrm{2}}^{\mathrm{\prime \prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathit{U}}_{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\\ & \mathrm{=}& \mathrm{-}\mathrm{3}\mathrm{+}\frac{{\mathit{C}}_{\mathrm{2}}}{\mathrm{2}}{\left(\frac{\mathrm{5}}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{-}\mathrm{1}\mathrm{+}\mathrm{3}{\cos }^{\mathrm{2}}{\mathit{\theta }}^{\mathrm{)}}\mathit{,}\end{array}$(37)where the quantity ${\mathit{C}}_{\mathrm{2}}\mathrm{\equiv }\frac{\mathrm{18}}{\mathrm{5}\mathit{\rho ̂}\mathrm{0}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{r̂}\mathrm{\prime }\frac{\mathrm{1}}{\mathit{r̂}\mathrm{\prime }}\mathit{\rho ̂}\mathrm{2}\mathrm{\left(}\mathit{r̂}\mathrm{\prime }\mathrm{\right)}$(38)depends implicitly on the parameter Ψ through the function and . The quantity C₂ is negative-definite since the quadrupole radial function of the density is negative on the entire domain of definition of the solution of the Poisson equation; the sign of the quadrupolar function is negative because the configurations in our family of models are always oblate. In passing, we also note that the the expansion around of the escape energy up to second order in radius is given by $\mathit{\psi }\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\mathrm{\Psi }\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }\mathit{r̂}\mathrm{2}}|}_{\mathrm{0}}\mathit{r̂}\mathrm{2}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{4}\mathrm{\right)}\mathit{.}$(39)Since the boundary of a configuration is defined by the condition , the density profile in the outer parts can be evaluated by performing an expansion with respect to ψ ≪ 1 $\mathit{\rho ̂}\mathrm{WT}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{=}\frac{\mathrm{9}}{\mathrm{5}}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }\hspace{0.17em}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathrm{;}$(40)the terms with fractional powers of ψ that appear in Eq. (34) cancel out with the first terms of the expansion of the double integral (which reduces to the incomplete gamma function).

The density profiles evaluated on the principal axes for a sequence of models with Ψ = 2, , and increasing values of the rotation strength parameter χ are illustrated in Fig. 11; the corresponding dimensionless escape energy profiles are displayed in Fig. 12. Configurations characterized by moderate rotation show monotonically decreasing profiles, whereas models with rapid rotation have the maximum value of the density profile in a position displaced with respect to the origin; in the extreme rotation regime, also the maximum value of the escape energy is off-centered. The sections of the isodensity and equipotential surfaces (presented in the first two rows of Fig. 13) clearly show that the offset of the density peak corresponds to the existence of a curious toroidal structure; the condition for the existence of such structure is discussed in Sect. 4.3.

4.2. The intrinsic kinematics

The calculation of the first order moment in velocity space of the distribution function ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ confirms that only the azimuthal component of the mean velocity is nonvanishing. The mean velocity in dimensionless form $\begin{array}{ccc}{\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}}_{\mathrm{WT}}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{3}}{{\mathrm{2}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{\rho ̂}\mathrm{WT}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\psi }}\mathrm{d}\mathit{s}\hspace{0.17em}\mathit{s}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathrm{d}\mathit{tt}\mathrm{\left[}\mathit{g}\mathrm{\right(}\mathit{s,t,}\mathit{r̂}\mathit{,\theta }\mathrm{)}{\mathrm{e}}^{\mathrm{-}\mathit{s}}{\mathrm{e}}^{\mathit{\psi }}\\ & & \mathrm{-}\ln \mathit{g}\mathrm{\left(}\mathit{s,t,}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{]}\end{array}$(41)can be calculated by numerically. As for the density profile, it is useful to evaluate the asymptotic expansion of Eq. (41) in the central regions and in the outer part of a given configuration. By performing a first order expansion in the radius with respect to the origin, we found that the mean rotation velocity reduces to the following expression ${\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}}_{\mathrm{WT}}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\mathrm{3}{\mathit{\chi }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{r̂}\sin \mathit{\theta }\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{3}\mathrm{\right)}\mathit{,}$(42)which corresponds to rigid rotation, with dimensionless angular velocity ; the expression does not depend on the concentration parameter, consistent with the asymptotic properties of the integral I(E,J_z).

In the outer parts of the models the mean rotation velocity is of order and thus vanishes at the boundary. The relevant numerical coefficients depend on the value of the parameter c; for example, for c = 1 we have ${\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}}_{\mathrm{WT}}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathrm{(}\mathrm{2}\mathrm{+}\mathrm{6}\mathit{b̅}\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{+}\mathrm{9}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{)}}{\mathrm{21}\mathrm{\left(}{\mathit{\chi }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{r̂}\sin \mathit{\theta }\mathrm{\right)}}\hspace{0.17em}\mathit{\psi }\mathit{.}$(43)The mean rotation velocity profiles on the equatorial plane for a selected sequence of models are displayed in Fig. 14: as the value of the central rotation strength parameter increases, the mean velocity profile becomes steeper in the inner parts and the maximum value increases; the central part of the profiles is well approximated by rigid rotation, as in Eq. (42).

By evaluating the second-order moments in velocity space of the distribution function ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$, it can be easily shown that, in our coordinate system, the pressure tensor is diagonal and that the radial component is equal to the polar one . Therefore, in the following, only the nontrivial components are discussed. The radial and azimuthal components in dimensionless form are given by $\begin{array}{ccc}\mathit{p̂}\mathrm{W}\mathit{,}\hspace{0.17em}\mathit{rr}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{3}}{\mathrm{4}}{\mathrm{e}}^{\mathit{\psi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\psi }}\mathrm{d}\mathit{s}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{s}}{\mathit{s}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathrm{d}\mathit{t}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{t}}^{\mathrm{2}}\mathrm{)}\mathit{g}\mathrm{\left(}\mathit{s,t,}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\\ & & \mathrm{-}\frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{-}\frac{\mathrm{4}}{\mathrm{35}}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\mathit{,}\\ \mathit{p̂}\mathrm{W}\mathit{,}\hspace{0.17em}\mathit{\phi \phi }\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{3}}{\mathrm{2}}{\mathrm{e}}^{\mathit{\psi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\psi }}\mathrm{d}\mathit{s}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{s}}{\mathit{s}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathrm{d}\mathit{t}{\mathit{t}}^{\mathrm{2}}\mathit{g}\mathrm{\left(}\mathit{s,t,}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\\ & & \mathrm{-}\frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{-}\frac{\mathrm{4}}{\mathrm{35}}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\mathrm{-}{\mathit{\rho ̂}}_{\mathrm{WT}}\mathrm{⟨}\mathit{v̂}\mathit{\phi }{\mathrm{⟩}}_{\mathrm{WT}}^{\mathrm{2}}\mathit{,}\end{array}$where the presence of the terms with fractional powers of ψ should be interpreted as in Eq. (34). The expansion up to second order in radius gives $\begin{array}{ccc}\mathit{p̂}\mathrm{W}\mathit{,}\hspace{0.17em}\mathit{rr}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,}\mathrm{\Psi }\mathrm{\right)}& \mathrm{=}& \mathit{\rho ̂}\mathit{WT,}\mathrm{0}\mathrm{+}\frac{\mathrm{1}}{\mathrm{5}}{\mathrm{e}}^{\mathrm{\Psi }}\mathit{\gamma }\left(\frac{\mathrm{7}}{\mathrm{2}}\mathit{,}\mathrm{\Psi }\right)\mathrm{[}\mathrm{9}\mathit{\chi }{\sin }^{\mathrm{2}}{\mathit{\theta }}^{}\\ & & \mathrm{+}{\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }\mathit{r̂}\mathrm{2}}|}_{\mathrm{0}}]\mathit{r̂}\mathrm{2}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{4}\mathrm{\right)}\mathit{.}\end{array}$(46)We also found that, to second order in radius, the pressure is isotropic . The components of the pressure tensor evaluated at the origin reduces to , consistent with the value obtained in the nonrotating limit (see Eq. (A.5) in Appendix A).

The asymptotic behavior of the pressure tensor components in the outer part can be written as $\begin{array}{ccc}\mathit{p̂}\mathrm{W}\mathit{,}\hspace{0.17em}\mathit{rr}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{18}}{\mathrm{35}}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }\hspace{0.17em}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{9}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathit{,}\\ \mathit{p̂}\mathrm{W}\mathit{,}\hspace{0.17em}\mathit{\phi \phi }\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{54}}{\mathrm{35}}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }\hspace{0.17em}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{9}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathit{,}\end{array}$respectively.

The relation between the dimensionless pressure and velocity dispersion tensor is given by $\mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathit{ij}\end{array}\mathrm{=}\mathit{p̂}{}_{\mathit{ij}}\mathit{/}\mathit{\rho ̂}$. The profiles of the radial and azimuthal component of the velocity dispersion tensor of a selected sequence of models, evaluated on the equatorial plane, are displayed in Fig. 15. For configurations in the moderate rotation regime the profiles are monotonically decreasing with the radius, with some variations in the slope in the intermediate and external parts, while for configurations in the rapid and extreme rotation regimes, their peak is off-centered. The sections in the meridional plane of the isobaric surfaces, defined with respect to the trace of the pressure tensor, show the presence of a central toroidal structure for the fast rotating models of the sequence (see the last row in Fig. 13).

The intrinsic kinematics can be further characterized by means of the anisotropy parameter⁴, defined as $\mathit{\alpha }\mathrm{=}\mathrm{1}\mathrm{-}\frac{{\mathit{p}}_{\mathit{\phi \phi }}}{{\mathit{p}}_{\mathit{rr}}}\mathrm{=}\mathrm{1}\mathrm{-}\frac{{\mathit{\sigma }}_{\mathit{\phi \phi }}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{rr}}^{\mathrm{2}}}\mathrm{·}$(49)From the asymptotic behavior of the pressure tensor components in the central region (see Eq. (46)), we find that α → 0 for , while from the expansion in outer parts (see Eqs. (47), (48)), we find that α →  − 2 as we approach the boundary. These limiting values do not depend on the dimensionless parameters that characterize the family ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$. In other words, the central region of the configurations is always characterized by isotropy in velocity space, while the regions next to the boundary show a strong tangentially-biased pressure anisotropy. The radial profiles of the anisotropy parameter for a selected sequence of models are displayed in Fig. 16. The values of α in the intermediate region of a given configuration depend on the values of the relevant dimensionless parameters . In fact, from an exploration of the entire four-dimensional parameter space, we found that, by increasing the value of c, the portion of a model dominated by radial anisotropy becomes slightly less extended, while, if the value of is increased, it increases (see Fig. 17, top panels). In turn, by increasing the value of the concentration, as measured by Ψ, the models are characterized by a significant radially-biased pressure anisotropy, which appears also in the intermediate region. By decreasing the value of the concentration, the intermediate region turns out to be dominated by tangential anisotropy, like in the outer parts (Fig. 17, bottom panels).

4.3. The condition for the existence of the central toroidal structure

In general, configurations with rapid rotation may exhibit highly deformed morphologies. In the context of rotating fluids, the regime of strong differential rotation has been successfully explored, at least for polytropes. Stoeckly (1965) and Geroyannis (1990) found that rapidly rotating polytropic fluids show a central toroidal structure; in addition, a self-consistent method for the construction of rapidly differentially rotating fluid systems with a great variety of shapes has been proposed by Hachisu (1986). In stellar dynamics, this regime has been rarely explored: Lynden-Bell (1962) and Prendergast & Tomer (1970) noted that some models with strong differential rotation show the density peak in a ring on their plane of symmetry.

The existence of a central toroidal structure in a given model of our family can be studied by using the asymptotic expansion of the density in the central regions, expressed in Eq. (36). In fact, if the second order derivative of the density with respect to the radius is positive, then the maximum value of the radial density profile is displaced from the geometric center of the configuration, so that a central toroidal structure is formed. Therefore, the relevant condition on the sign of the derivative can be expressed as $\mathrm{9}\mathit{\chi }{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{-}\mathrm{3}\mathrm{+}\frac{{\mathit{C}}_{\mathrm{2}}}{\mathrm{2}}{\left(\frac{\mathrm{5}}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{-}\mathrm{1}\mathrm{+}\mathrm{3}{\cos }^{\mathrm{2}}{\mathit{\theta }}^{\mathrm{)}}\mathit{>}\mathrm{0}\mathit{,}$(50)which, on the equatorial plane (i.e., at θ = π/2), reduces to a simple condition for the central rotation strength parameter $\mathit{\chi }\mathit{>}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{+}\frac{{\mathit{C}}_{\mathrm{2}}}{\mathrm{18}}{\left(\frac{\mathrm{5}}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(51)The expression on the right-hand side depends implicitly on the dimensionless parameters Ψ, , and c through the quantity C₂, defined by Eq. (38). Since C₂ is negative-definite, the condition χ > 1/3 provides the upper limit to Eq. (51). Actually, the models characterized by values of χ immediately above the threshold given by Eq. (51), show a very small central toroidal structure, with a shallow increase of the density with respect to the geometric center of the configuration (e.g., for the model with Ψ = 2 and χ = 0.36, illustrated in Fig. 11, the density peak, at the center of the toroidal structure, is merely log (ρ/ρ₀) = 0.014 and the central structure itself is barely visible in the meridional sections of the isodensity surfaces, depicted in the third panel of the first row of Fig. 13). To construct a model with a sizable central toroidal structure, the value of the parameter χ should be at least twice the threshold given by the above-mentioned condition (e.g., see the last panel of the first row of Fig. 13). In conclusion, by comparing Figs. 9 and 18, we note that, for configurations in the moderate rotation regime (i.e., with ), the central toroidal structure is absent or very small, while starting from the rapid rotation regime, the structure is always present and becomes more extended as the value of the central rotation strength increases.

The presence of the central toroidal structure can be also be interpreted in terms of to the intrinsic kinematical properties the models. From the radial component of the Jeans equation, expressed in dimensionless spherical coordinates, $\frac{\mathrm{1}}{\mathit{\rho ̂}}\frac{\mathit{\partial }\mathit{\rho ̂}}{\mathit{\partial }\mathit{r̂}}\frac{\mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathit{rr}\end{array}}{\mathit{r̂}}\mathrm{=}\frac{\mathrm{1}}{\mathit{r̂}}\frac{\mathit{\partial \psi }}{\mathit{\partial }\mathit{r̂}}\mathrm{+}\frac{\mathrm{⟨}\mathit{v̂}\mathit{\phi }{\mathrm{⟩}}^{\mathrm{2}}}{\mathit{r̂}\mathrm{2}}\mathrm{-}\left[\frac{\mathrm{1}}{\mathit{r̂}}\frac{\mathit{\partial }\mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathit{rr}\end{array}}{\mathit{\partial }\mathit{r̂}}\mathrm{+}\frac{\mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathit{rr}\end{array}\mathrm{-}\mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathit{\phi \phi }\end{array}}{\mathit{r̂}\mathrm{2}}\right]\mathit{,}$(52)the sign of the first order derivative of the density with respect to the radius (on the left-hand side of Eq. (52)) depends on (i) the difference between the angular velocity associated with the circular orbit of a single star and the angular velocity of the model , associated to the mean rotation velocity , (ii) a more complex pressure term (in square brackets on the right-hand side of Eq. (52)). To check for the presence of the central toroidal structure, it is sufficient to study the Jeans equation in the central region of the model. Therefore, by inserting the relevant asymptotic expansions for the mean velocity and the escape energy (see Eqs. (42) and (39)), the term (i) reduces to the expression indicated on the left-hand side of Eq. (50), and, by inserting the relevant asymptotic expansions for the pressure tensor components, the term (ii) can be written as $\begin{array}{ccc}& & \left[\mathrm{9}\mathit{\chi }{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{-}\mathrm{3}\mathrm{+}\frac{{\mathit{C}}_{\mathrm{2}}}{\mathrm{2}}{\left(\frac{\mathrm{5}}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{-}\mathrm{1}\mathrm{+}\mathrm{3}{\cos }^{\mathrm{2}}{\mathit{\theta }}^{\mathrm{)}}\right]\\ & & \mathrm{\times }\hspace{0.17em}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{5}}{\mathrm{7}}\frac{\mathit{\gamma }\mathrm{(}\mathrm{9}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}\mathit{\gamma }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}}{\mathit{\gamma }\mathrm{(}\mathrm{7}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }{\mathrm{)}}^{\mathrm{2}}}\right]\mathrm{·}\end{array}$(53)By combining the asymptotic expressions of terms (i) and (ii), it follows⁵ that the requirement of a positive density gradient on the left-hand side of the Jeans equation is equivalent to the condition expressed by Eq. (50).

In conclusion, we have independently tested the validity of the condition for the existence of the central toroidal structure derived at the beginning of this section. We have also shown that the requirement of positivity of term (i) in the radial Jeans equation is a necessary and sufficient condition for the presence of the central toroidal structure. In other words, the central toroidal structure exists if and only if, in the central region, the angular velocity associated with the internal rotation is higher than the angular velocity associated with the circular orbit of a single star. Figure 19 shows the relevant angular velocities for the first three models of the sequence considered in Fig. 13; it is apparent that, for the configuration in which the central toroidal structure is present, . This result strictly depends on the adopted truncation in phase space for the distribution function ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$; in fact, in Appendix B we show that for the alternative distribution function ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ , the condition on the angular velocities is a necessary but not sufficient condition for the existence of the central toroidal structure.

4.4. The condition for maximally rotating models

If we consider a sequence of models with fixed values of and increasing values of χ, the central toroidal structure becomes progressively more extended and characterized by a larger aperture angle (i.e., the angle spanned at small radii by the isodensity surface that represents the boundary of the central toroidal structure). For high values of the central rotation strength parameter χ, a smaller central toroidal structure appears also in the equipotential surfaces (see the escape energy profile of the fastest rotating model in Fig. 12). These two structures can be characterized in terms of θ_p and θ_d, defined as the complements of the semi-aperture angle of the inner cusp of toroidal structures in the equipotential and isodensity surface, respectively. In fact, since the boundary of the central cusp of the toroidal structures, is defined by , and , by using the relevant asymptotic expansion up to second order in given in Eqs. (39) and (36), the following expressions for the angles are obtained: $\begin{array}{ccc}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{p}}& \mathrm{=}& \frac{\mathrm{3}\mathrm{+}\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{C}}_{\mathrm{2}}\mathit{/}\mathrm{2}}{\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{C}}_{\mathrm{2}}}\mathit{,}\\ {\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{d}}& \mathrm{=}& \frac{\mathrm{3}\mathrm{+}\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{C}}_{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{-}\mathrm{9}\mathit{\chi }}{\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{C}}_{\mathrm{2}}\mathrm{-}\mathrm{9}\mathit{\chi }}\mathrm{·}\end{array}$The two angles are related, because ${\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{d}}\mathrm{=}\frac{\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{C}}_{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{p}}\mathrm{-}\mathrm{9}\mathit{\chi }}{\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{C}}_{\mathrm{2}}\mathrm{-}\mathrm{9}\mathit{\chi }}\mathrm{·}$(56)Note that θ_p and θ_d decrease as the values of  | C₂ |  and χ increase (we recall that C₂ is negative-definite), that is the corresponding semi-aperture angles become larger.

For high values of the central rotation strength χ, which imply a high degree of quadrupolar deformation, as measured by the quantity  | C₂ | , it is easy to see that cos²θ_p → 1/3; interestingly, we also found (numerically) that a limiting value exists also for θ_d, given by cos²θ_d → 2/3. By inserting the limiting value of θ_p in Eq. (56), the condition for the maximum value of the angle θ_d (i.e., cos²θ_d < 2/3) can be translated in a simple condition for the central rotation strength parameter $\mathit{\chi }\mathit{<}\frac{\mathrm{1}}{\mathrm{6}}{\left(\frac{\mathrm{5}}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\left|}{\mathit{C}}_{\mathrm{2}}\mathrm{\right|}\mathit{,}$(57)which basically provides a limit to the deformation induced by the central rotation itself. The previous condition has been checked numerically and we found that, for values of χ above the threshold, the iteration scheme used for the solution of the Poisson equation does not converge.

5. Projected properties of the differentially rotating models

5.1. The surface density profile

The calculation of the projected properties has been performed by following the same projection rules described in Sect. 2.4, that is the line of sight corresponds to the ẑ_P-axis of a new frame of reference in which the projection plane is denoted by . In particular, we studied in detail the “edge-on” view, that is the projection along the -axis of the intrinsic frame of reference (, , and ŷ_P = ẑ, i.e., the viewing angles of the rotation matrix are θ = π/2, φ = 0). The dimensionless surface density distribution has been calculated by numerical integration of Eq. (19) on an equally-spaced square cartesian grid on the projection plane.

The first row of Fig. 20 shows the contour maps of the surface density of selected differentially rotating models in the moderate and rapid rotation regime. As result of projection, the dimples on the rotation axis, which are prominent in the corresponding intrinsic isodensity surfaces (see Fig. 13 for the meridional sections of the same sequence of models), are less pronounced. In addition, the central toroidal structure in the projected density distribution is visible only if it has a reasonable size in the intrinsic density distribution (see third and last model of the sequence illustrated in Fig. 20).

The surface density profiles of the same sequence of differentially rotating models, evaluated along the principal axes of the projection plane, are presented in Fig. 21. For the configurations in which the central toroidal structure in projection is absent, we calculated the relevant ellipticity profiles, as functions of the semi-major axis of the projected image â_P, and they are reported in Fig. 22. As expected, the configurations in the moderate rotation regime are characterized by nonmonotonic ellipticity profiles, while models in the rapid rotation regime have monotonically decreasing profiles; to some extent, this morphological feature is complementary to that of the uniformly rotating models, in which the configurations are always characterized by monotonically increasing ellipticity profiles. Interestingly, the behavior of the ellipticity profiles does not necessarily correlate with the mean line-of-sight velocity profiles (see next subsection), that is configurations with a nonmonotonic velocity profile may have a monotonic ellipticity profile.

In addition, the isophotes of models show clear departures from a pure ellipse, which, at variance with the family of rigidly rotating models, can be characterized as a “boxy” overall trend, particularly evident in the intermediate parts of the configurations (see third and fourth panels of the first row of Fig. 20).

5.2. The line-of-sight velocity distribution

The projected velocity moments are calculated by integrating along the line of sight (weighted by the intrinsic density) the corresponding intrinsic quantities. As for the surface density distribution, we studied in detail the kinematics projected along the -axis of the intrinsic coordinate system. Therefore, the dimensionless mean line-of-sight velocity and the line-of-sight velocity dispersion can be written as $\begin{array}{ccc}& & \mathrm{⟨}\mathit{v̂}\mathrm{P}\mathrm{⟩}\mathrm{\left(}\mathit{x̂}\mathrm{P}\mathit{,}\mathit{ŷ}\mathrm{P}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{\Sigma ̂}\mathrm{\left(}\mathit{x̂}\mathrm{P}\mathit{,}\mathit{ŷ}\mathrm{P}\mathrm{\right)}}{\mathrm{\int }}_{\mathit{ẑ}\mathrm{sp}}^{\mathit{ẑ}\mathrm{sp}}\mathrm{d}\mathit{ẑ}\mathrm{P}\hspace{0.17em}\mathit{\rho ̂}\mathrm{\left(}\hat{{r}}\mathrm{P}\mathrm{\right)}\hspace{0.17em}\frac{\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}\hspace{0.17em}\mathit{x̂}\mathrm{P}}{\mathrm{(}\mathit{ẑ}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}\mathrm{+}\mathit{x̂}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{,} \\ & & \mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}\mathrm{\left(}\mathit{x̂}\mathrm{P}\mathit{,}\mathit{ŷ}\mathrm{P}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Sigma ̂}\mathrm{\left(}\mathit{x̂}\mathrm{P}\mathit{,}\mathit{ŷ}\mathrm{P}\mathrm{\right)}}{\mathrm{\int }}_{\mathit{ẑ}\mathrm{sp}}^{\mathit{ẑ}\mathrm{sp}}\mathrm{d}\mathit{ẑ}\mathrm{P}\hspace{0.17em}\mathit{\rho ̂}\mathrm{\left(}\hat{{r}}\mathrm{P}\mathrm{\right)}\hspace{0.17em}[\frac{\mathrm{1}}{\mathit{ẑ}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}\mathrm{+}\mathit{x̂}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}}\mathrm{(}\mathit{ẑ}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}\mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathit{rr}\end{array}\\ & & \mathrm{+}\mathit{x̂}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}\mathit{\sigma ̂}\begin{array}{}\mathrm{2}\\ \mathit{\phi \phi }\end{array}\mathrm{+}\mathrm{⟨}\mathit{v̂}\mathit{\phi }{\mathrm{⟩}}^{\mathrm{2}}\mathit{x̂}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}\mathrm{)}\mathrm{+}\mathrm{⟨}\mathit{v̂}\mathrm{P}{\mathrm{⟩}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{⟨}\mathit{v̂}\mathrm{P}\mathrm{⟩}\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}\frac{\mathit{x̂}\mathrm{P}}{\mathrm{(}\mathit{ẑ}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}\mathrm{+}\mathit{x̂}\begin{array}{}\mathrm{2}\\ \mathrm{P}\end{array}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}]\mathrm{·} \end{array}$Selected mean line-of-sight velocity profiles (for the sequence of models presented in Fig. 21), evaluated along the -axis of the projection plane, are illustrated in Fig. 23, where the sign of the mean velocity in the two half-planes is consistent with Eq. (58). As the value of the central rotation strength parameter increases, the slope of the inner part of the profile becomes steeper, since the asymptotic behavior of the intrinsic velocity in the central regions is approximately that of a rigid rotation, with the angular velocity proportional to χ^1/2 (see Eq.  (42)); in addition, because the entire configuration becomes more compact, the radial position of the peak of the velocity profile shrinks progressively.

The line-of-sight velocity dispersion profiles of the same sequence of models, evaluated along the principal axes of the projection plane, are presented in Fig. 24. For configurations in the moderate rotation regime, the variations in the slope at intermediate and outer radii, which characterize the azimuthal component of the intrinsic velocity dispersion tensor, are still visible in projection, while the inner part of the profile has a flat core. For models in the rapid and extreme rotation regime, the maximum value of the profile is displaced from the geometric center; in this case, the inner part of the profile has a nontrivial gradient.

We also calculated the ratio of the mean line-of-sight (rotation) velocity to the central line-of-sight velocity dispersion, as a measure of the amount of ordered motions compared to random motions (see Fig. 25 for relevant profiles evaluated along the -axis, for the sequence of models discussed above). The models in the moderate rotation regime show values that are consistent with those observed in Galactic globular clusters (see Introduction). The configurations in the rapid and extreme rotation regime are characterized by higher values of the ratio, that can be even greater than one; such high values are measured only in the class of elliptical galaxies known as fast rotators (see Davies et al. 1983). Of course, the values of this ratio strongly depend on the line of sight on which the projection is performed (we recall that here we illustrate the results obtained only from the “edge-on” view, which is the most favorable for the detection of ordered motions).

6. Effect of truncation in phase space

To explore the nontrivial effects of the two options for the truncation prescription of the family of differentially rotating models (see Sect. 3.1; see also Hunter 1977), we selected two representative models and studied the relevant phase space density N(Ê,Ĵ_z), defined in such a way that the total dimensionless mass is given by , with the specific energy Ê = a   E and z-component of the specific angular momentum $\widehat{{\mathit{J}}_{\mathit{z}}}\mathrm{=}\mathrm{(}\hat{{r}}\mathrm{\times }{v̂}{\mathrm{)}}_{\mathit{z}}\mathrm{=}\mathit{v̂}{}_{\mathit{\phi }}\mathit{r̂}\sin \mathit{\theta }$. The histograms of the phase space density, constructed by means of a Monte Carlo sampling of the distribution functions, are presented in the top panels of Fig. 26.

We also constructed the corresponding Lindblad diagrams (Lindblad 1933), that is the representation in the plane (Ê,Ĵ_z) of orbits in phase space in a given model (see bottom panels of Fig. 26). For any given Ĵ_z, the orbit with the lowest energy corresponds to a circular orbit in the equatorial plane; therefore, in both families of models, the circular orbits form the lower cuspy boundary of the model in the diagram. The upper boundary corresponds to the energy truncation, introduced by the cut-off constant Ê₀. Of course, since both families of models are characterized by internal rotation, the distribution of the orbits is asymmetric with respect to the z-component of the angular momentum, but different truncation prescriptions determine a different distribution of the orbits in phase space. In particular, the Wilson truncation introduces a sharp depopulation of retrograde orbits with intermediate energy, while the plain truncation is associated with a smooth decrease of the number of retrograde orbits with high energy. Alternative options for the truncation prescription in rotating models correspond to a distribution of orbits in different regions of the diagram (see Fig. 1 in Rowley 1988). Note that a truncation with respect to the Jacobi integral H = E − ωJ_z only, as in our family of rigidly rotating models ${\mathit{f}}_{\mathrm{K}}^{\mathit{r}}\mathrm{\left(}\mathit{H}\mathrm{\right)}$ (see Eq. (1)), corresponds to a region in the diagram bounded by a straight line (with ω as slope) and the curves corresponding to the energy of circular orbits.

In addition, as the concentration of the models increases, the cusp of the lower boundary, given by the energy of circular orbits, becomes sharper; as noted also by Rowley (1988), this is the reason why centrally concentrated configurations with rigid rotation cannot be very flat, except for the outer parts, as described in Sect. 2.4.

7. Discussion and conclusions

In this paper we have constructed two new families of self-consistent axisymmetric models of quasi-relaxed stellar systems, characterized by the presence of internal rotation (see Table 1); a full description in terms of the intrinsic and projected properties has been provided. The main results can be summarized as follows:

• Driven by general statistical mechanics considerations, we started by constructing a family of rigidly rotating dynamical models; this family is defined as an extension of the King (1966) models to the case of axisymmetric equilibria flattened by solid-body rotation, with the relevant distribution function dependent only on the Jacobi integral. The configurations have been constructed self-consistently by solving the Poisson-Laplace equation for the mean-field potential by means of a perturbation method (described in Paper I), using a measure of the rotation strength as the expansion parameter. The two-dimensional parameter space which characterizes the family (concentration and rotation strength) can be described in terms of two regimes. Models in the low-deformation regime are almost indistinguishable from the corresponding spherical King models. Highly-deformed models are quasi-spherical in the central regions and show significant deviations from spherical symmetry in the outer parts; in particular, they are flattened toward the equatorial plane and exhibit a sort of “disky” appearance. The resulting eccentricity profile is a monotonically increasing function of radius; the (finite) central value can be expressed analytically in terms of the rotation strength parameter. From the kinematical point of view, the models are characterized by pressure isotropy and cylindrical rotation.

  Fig. 26 Top panels: histogram (scaled to unity) of the phase space density N(Ê,Ĵz) for a differentially rotating model characterized by Wilson truncation (left) with , and for one characterized by plain truncation (right) with . Both models are characterized by . Bottom panels: Lindblad diagrams for the same models presented in the top panels. The graphs have been obtained by a Monte Carlo sampling of the relevant distribution functions with 65   536 particles.

• In view of possible applications to globular clusters, we have constructed a second family of dynamical models, characterized by differential rotation, designed to be approximately rigid in the central regions and to vanish in the outer parts, where the imposed energy truncation is effective. In this case, the relevant Poisson equation is solved by means of a spectral iteration method based on the Legendre expansion of the density and the potential. The full parameter space is now four-dimensional, with two additional parameters, defining the shape of the rotation profile. Three rotation regimes can be introduced, namely of moderate, rapid, and extreme rotation. However, significant variations in the structure of the models are primarily associated with concentration and central rotation strength, as for the previous family. We explored the properties of the configurations resulting from two options for the truncation prescription, with emphasis on the family which, in the limit of vanishing internal rotation, reduces to the spherical limit of the models proposed by Wilson (1975). In particular, configurations in the rapid and extreme rotation regimes exhibit a central toroidal structure, the volume of which increases with the value of the central rotation strength parameter. By making use of asymptotic expansions of the density, mean velocity, and pressure tensor components for small radii, we found the condition for the existence of such central toroidal structure, as well as the condition for the maximum value of the central rotation strength parameter admitted by a configuration with a given concentration.

• The differentially rotating models show a variety of realistic velocity dispersion profiles, characterized by the presence of pressure isotropy and radially-biased anisotropy in the central and intermediate regions, respectively. The kinematical behavior in the outer parts depends on the adopted truncation prescription; in particular, the family which, in the nonrotating limit, reduces to the Wilson spheres is characterized by tangentially-biased anisotropy. This kinematical feature (rarely obtained in equilibrium models) is of great interest for two reasons: (i) tangentially-biased pressure anisotropy is observed in the presence of internal rotation in globular clusters. For example, the full three-dimensional view of the velocity space of ω Cen, obtained from proper motions and radial velocities measurements, has revealed that this object is characterized by significant rotation and tangential anisotropy in the outer parts (van de Ven et al. 2006); (ii) the dynamical evolution of a cluster in a tidal field is known to induce a rapid development of tangential anisotropy in the outer parts of the stellar system. In fact, if a cluster fills its Roche lobe and starts losing mass, there is a preferential loss of stars on radial orbits induced by the external tidal field at large radii, where tangential anisotropy in velocity space is thus established (Takahashi & Lee 2000; Baumgardt & Makino 2003).

• The presence of differential rotation may induce nontrivial gradients in the line-of-sight velocity dispersion profile of a stellar system, even if the amount of rotation is modest. Therefore, this important physical ingredient should be taken into account properly. In this respect, dynamical studies of globular clusters and other low-mass stellar systems by means of models based on the use of the Jeans equations are less satisfactory, because, at least in their most popular (nonrotating) application, they are used to reproduce variations in the slope of the kinematical profile of a system only by means of a (sometimes significant) amount of pressure anisotropy.

• As expected, differential rotation also induces nontrivial deviations from spherical symmetry; in fact, the models are characterized by a great variety of (projected) ellipticity profiles, dependent on the combined effect of concentration and central rotation strength. Configurations in the moderate rotation regime are characterized by realistic nonmonotonic ellipticity profiles (e.g., see Geyer et al. 1983), while models in the rapid rotation regime have monotonically decreasing profiles. To some extent, this morphological feature is complementary to that of the uniformly rotating models, in which the configurations are always characterized by monotonically increasing ellipticity profiles. Interestingly, the behavior of the ellipticity profiles is not necessarily correlated with the mean line-of-sight velocity profiles, that is configurations with a nonmonotonic mean velocity profile may have a monotonic ellipticity profile. In addition, the isophotes of the relevant surface density distribution tend to be characterized by a “boxy” structure.

• From a comparison of the equilibrium configurations resulting from two options for the truncation prescription of the family of differentially rotating models, we confirm that the interplay between internal rotation, anisotropy in velocity space, and truncation in phase space is highly nontrivial. In fact, as also noted by Hunter (1977), the structure of the outer parts of a model is particularly sensitive, both from the morphological and the kinematical point of view, to the adopted truncation. One way to select the most appropriate truncation from the physical point of view will be to address the issue in the context of formation and evolution of the class of stellar systems under consideration (see also last item below).

• Models in the moderate rotation regime seem to be particularly appropriate for describing rotating globular clusters, since the relevant configurations are characterized by a number of realistic properties, such as the presence of nonmonotonic ellipticity profile, the behavior of surface density profile in the outer parts similar to the one associated with spherical Wilson models, the existence of pressure isotropy in the central regions and tangentially-biased anisotropy at the boundary, as well as realistic values of the ratio  ⟨ v_P ⟩ /σ_P,   0. In a following paper, we plan to apply our family of differentially rotating models to selected Galactic globular clusters that show the presence of significant rotation, such as ω Cen and 47 Tuc.

• Configurations with strong differential rotation, characterized by the presence of a sizable central toroidal structure and by a off-centered peak of the surface brightness profiles, may be useful to shed light on the internal dynamics of the so-called “ring clusters”. This class of object, originally observed in the Small Magellanic Cloud (Hill & Zaritsky 2006) and subsequently noted also in the Large Magellanic Cloud (Werchan & Zaritsky 2011), is characterized by a sizable dimple of the central surface brightness, resulting in an off-centered peaked density profile. A proper dynamical interpretation of these objects is currently missing.

• The families of models illustrated in the present paper may also help to clarify the role of angular momentum in the formation and dynamical evolution of globular clusters. The results of an extensive survey of N-body simulations, designed to study the dynamical stability and the long term evolution of the models described here, will be presented in follow-up papers (Varri et al. 2011, and in prep.).

---

Acknowledgments

We would like to thank E. Vesperini for a critical reading of the manuscript and a number of valuable suggestions, and S. L. W. McMillan, D. Heggie, L. Ciotti, J. Anderson, F. Rasio, G. Lodato, A. Sollima, P. Bianchini, and A. Zocchi for interesting conversations. A.L.V. gratefully acknowledges the hospitality of the Department of Physics at Drexel University, where the first part of the manuscript was written, under the support of a US-Italy Fulbright Visiting Student Researcher Fellowship. This work is partly supported by the Italian MIUR.

References

Appendix A: Spherical limit of rotating models

The spherical nonrotating limit of the families of rotating models considered in the present paper are defined by the following distribution functions $\begin{array}{ccc}& & {\mathit{f}}_{\mathrm{PT}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\mathit{A}{\mathrm{e}}^{\mathrm{-}\mathit{a}{\mathit{E}}_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{-}\mathit{a}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}_{\mathrm{0}}\mathrm{)}}\mathit{,}\\ & & {\mathit{f}}_{\mathrm{K}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\mathit{A}{\mathrm{e}}^{\mathrm{-}\mathit{a}{\mathit{E}}_{\mathrm{0}}}\mathrm{[}{\mathrm{e}}^{\mathrm{-}\mathit{a}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}_{\mathrm{0}}\mathrm{)}}\mathrm{-}{\mathrm{1}}^{\mathrm{]}}\mathit{,}\\ & & {\mathit{f}}_{\mathrm{WT}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\mathit{A}{\mathrm{e}}^{\mathrm{-}\mathit{a}{\mathit{E}}_{\mathrm{0}}}\mathrm{[}{\mathrm{e}}^{\mathrm{-}\mathit{a}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}_{\mathrm{0}}\mathrm{)}}\mathrm{-}\mathrm{1}\mathrm{+}\mathit{a}\mathrm{(}\mathit{E}\mathrm{-}{\mathit{E}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{]}}\mathit{,}\end{array}$if E ≤ E₀ and f_i(E) = 0 otherwise (in the following, the index i = 1,2,3 denotes the models defined by Eqs. (A.1)–(A.3), in the same order). In particular, f_K(E) defines the King (1966) models and represents the spherical nonrotating limit of the family of rigidly rotating models defined by ${\mathit{f}}_{\mathrm{K}}^{\mathrm{r}}\mathrm{\left(}\mathit{H}\mathrm{\right)}$ (see Eq. (1)); f_WT(E) and f_PT(E) are the spherical isotropic limit of the Wilson (1975) and Prendergast & Tomer (1970) models, and represent the spherical nonrotating limit of our differentially rotating models defined by ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ (see Eq. (21)) and ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ (see Eq. (22)), respectively. In the previous expressions, A, a, are positive constants, defining two dimensional scales, while E₀ is the cut-off energy, which implies the existence of a truncation radius r_tr for the spherical system. For all the families of models considered here one important parameter is the central concentration, as measured by the depth of the dimensionless central potential well Ψ = ψ(0) = a [E₀ − Φ(0)]  or by the parameter . Note that, if we compare three models (one for each family) having the same value of Ψ, the corresponding values of c are not the same, that is the relevant dimensionless truncation radii are different ( increases from f_PT(E), to f_K(E), to f_WT(E)). In other words, the structure of the outer parts of a spherical isotropic truncated model strictly depends on the truncation prescription; this property is particularly relevant for the interpretation of the photometric profiles of globular clusters (see the systematic comparison between spherical King and Wilson models performed by McLaughlin & van der Marel 2005).

From the integration of the distribution functions in velocity space, the corresponding intrinsic density distributions are recovered. Using the same dimensionless units introduced in the main text, we denote the dimensionless density profiles by (for the definition of Â, see Eq. (9)), with $\mathit{\rho ̂}\mathrm{i}\mathit{,}\mathrm{S}\mathrm{\left(}\mathit{\psi }\mathrm{\right)}\mathrm{=}{\mathit{D}}_{\mathrm{i}}\hspace{0.17em}{\mathrm{e}}^{\mathit{\psi }}\mathit{\gamma }\left({\mathit{E}}_{\mathrm{i}}\mathit{,\psi }\right)\mathit{,}$(A.4)here ψ indicates the dimensionless escape energy, defined as in Eq. (24).

Similarly, the trace of the pressure tensor in dimensionless form (divided by a factor 3) can be written as , where $\mathit{p̂}\mathrm{i}\mathit{,}\mathrm{S}\mathrm{\left(}\mathit{\psi }\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathit{i}}\hspace{0.17em}{\mathrm{e}}^{\mathit{\psi }}\mathit{\gamma }\left({\mathit{E}}_{\mathit{i}}\mathrm{+}\mathrm{1}\mathit{,\psi }\right)\mathrm{;}$(A.5)D_i, P_i, and E_i are numerical coefficients resulting from the integration in velocity space and are summarized in Table A.1. Note that the coefficient appearing as first argument of the incomplete gamma function in the pressure profile is related to the corresponding coefficient in the density profile, because they are, respectively, the second-order and zeroth-order moment of the distribution function in velocity space.

Appendix B: Plain-truncated rotating models

We summarize here the intrinsic properties of the family of models defined by ${\mathit{f}}_{\mathrm{PT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ (see Eq. (22)), for a comparison with the corresponding quantities derived from the family of models defined by ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ (see Eq. (21)); a full description is given by Varri (2012). The density profile in dimensionless form can be written as $\mathit{\rho ̂}\mathrm{PT}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{=}\mathit{\rho ̂}\mathrm{WT}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{+}{\mathit{\psi }}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{;}$(B.1)the asymptotic behavior of the density profile in the outer parts (ψ → 0), with respect to the dimensionless escape energy, is given by: $\mathit{\rho ̂}\mathrm{PT}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{=}{\mathit{\psi }}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{9}}{\mathrm{2}}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }\right)\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathit{.}$(B.2)In the central region (), the expansion to second order in radius is $\begin{array}{ccc}\mathit{\rho ̂}\mathrm{PT}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,}\mathrm{\Psi }\mathrm{\right)}& \mathrm{=}& \mathit{\rho ̂}\mathrm{PT}\mathit{,}\mathrm{0}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}[\mathrm{9}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }{\mathrm{e}}^{\mathrm{\Psi }}\mathit{\gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}\mathit{,}\mathrm{\Psi }\right)\\ & & {\mathrm{+}\frac{\mathrm{3}}{\mathrm{4}}{\mathrm{e}}^{\mathrm{\Psi }}\mathit{\gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\mathrm{\Psi }\right)\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }\mathit{r̂}\mathrm{2}}|}_{\mathrm{0}}]\mathit{r̂}\mathrm{2}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{4}\mathrm{\right)}\mathit{,}\end{array}$(B.3)where the central value is , consistent with the value of the corresponding nonrotating model.

The mean velocity is in the azimuthal direction. In dimensionless units it is given by ${\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}}_{\mathrm{PT}}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}\hspace{0.17em}{\mathrm{e}}^{\mathit{\psi }}}{{\mathrm{2}}^{\mathrm{3}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathit{\rho ̂}\mathrm{PT}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\psi }}\mathrm{d}\mathit{s}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{s}}\mathit{s}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathrm{d}\mathit{t}\hspace{0.17em}\mathit{t}\hspace{0.17em}\mathit{g}\mathrm{\left(}\mathit{s,t,}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathit{,}$(B.4)which, in the outer parts of the models reduces to ${\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}}_{\mathrm{PT}}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{=}\mathit{\psi }\frac{\mathrm{6}}{\mathrm{5}}{\mathit{\chi }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{r̂}\sin \mathit{\theta }\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathit{.}$(B.5)The expansion for to first order in radius can be written as ${\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}}_{\mathrm{PT}}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,}\mathrm{\Psi }\mathrm{\right)}\mathrm{=}\mathrm{2}\frac{\mathit{\gamma }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}}{\mathit{\gamma }\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}}{\mathit{\chi }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{r̂}\sin \mathit{\theta }\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{3}\mathrm{\right)}\mathrm{;}$(B.6)at variance with the family defined by ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$, the dimensionless angular velocity depends also on the concentration parameter.

As far as the pressure tensor is concerned, we recall that by construction . The dimensionless radial component is given by $\mathit{p̂}\mathrm{PT}\mathit{,}\hspace{0.17em}\mathit{rr}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{=}\mathit{p̂}\mathrm{W}\mathit{,}\hspace{0.17em}\mathit{rr}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{+}\frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{+}\frac{\mathrm{4}}{\mathrm{35}}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\mathit{,}$(B.7)which, expanded in ψ (as is appropriate for the outer parts), reduces to $\begin{array}{ccc}\mathit{p̂}\mathit{PT,rr}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\\ & & \mathrm{+}\frac{\mathrm{4}}{\mathrm{35}}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{9}}{\mathrm{2}}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }\right)\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{9}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathit{.}\end{array}$(B.8)In the inner parts it can be approximated to second order in radius, by the following expression $\begin{array}{ccc}\mathit{p̂}\mathrm{PT}\mathit{,rr}\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,}\mathrm{\Psi }\mathrm{\right)}& \mathrm{=}& \mathit{p̂}\mathrm{PT}\mathit{,}\hspace{0.17em}\mathrm{0}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}[\frac{\mathrm{18}}{\mathrm{5}}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }{\mathrm{e}}^{\mathrm{\Psi }}\mathit{\gamma }\left(\frac{\mathrm{7}}{\mathrm{2}}\mathit{,}\mathrm{\Psi }\right)\\ & & {\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{\Psi }}\mathit{\gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}\mathit{,}\mathrm{\Psi }\right)\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }\mathit{r̂}\mathrm{2}}|}_{\mathrm{0}}]\mathit{r̂}\mathrm{2}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{4}\mathrm{\right)}\mathit{,}\end{array}$(B.9)where the central value is given by , consistent with the value of the corresponding model in the limit of vanishing rotation.

The azimuthal component of the pressure tensor can be written as $\begin{array}{ccc}\mathit{p̂}\mathrm{PT}\mathit{,}\hspace{0.17em}\mathit{\phi \phi }\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \mathit{p̂}\mathrm{W}\mathit{,}\hspace{0.17em}\mathit{\phi \phi }\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}\mathrm{+}\frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{+}\frac{\mathrm{4}}{\mathrm{35}}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\\ & & \mathrm{+}\mathit{\rho ̂}\mathrm{WT}{\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}}_{\mathrm{WT}}^{\mathrm{2}}\mathrm{-}\mathit{\rho ̂}\mathrm{PT}{\mathrm{⟨}\mathit{v̂}\mathit{\phi }\mathrm{⟩}}_{\mathrm{PT}}^{\mathrm{2}}\mathit{,}\end{array}$(B.10)which at the boundary is approximated by $\begin{array}{ccc}\mathit{p̂}\mathrm{PT}\mathit{,\phi \phi }\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,\psi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{2}}{\mathrm{5}}{\mathit{\psi }}^{\mathrm{5}\mathit{/}\mathrm{2}}\\ & & \mathrm{+}\frac{\mathrm{4}}{\mathrm{35}}{\mathit{\psi }}^{\mathrm{7}\mathit{/}\mathrm{2}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{9}}{\mathrm{10}}\mathit{\chi }\mathit{r̂}\mathrm{2}{\sin }^{\mathrm{2}}\mathit{\theta }\right)\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{\psi }}^{\mathrm{9}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathit{.}\end{array}$(B.11)Close to the center we find $\begin{array}{ccc}& & \mathit{p̂}\mathrm{PT}\mathit{,\phi \phi }\mathrm{\left(}\mathit{r̂}\mathit{,\theta ,}\mathrm{\Psi }\mathrm{\right)}\mathrm{=}\mathit{p̂}\mathrm{PT}\mathit{,}\mathrm{0}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\{\mathrm{6}\mathit{\chi }{\sin }^{\mathrm{2}}\mathit{\theta }[\frac{\mathrm{9}}{\mathrm{5}}{\mathrm{e}}^{\mathrm{\Psi }}\mathit{\gamma }\left(\frac{\mathrm{7}}{\mathrm{2}}\mathit{,}\mathrm{\Psi }\right)\\ & & \mathrm{-}{\mathrm{2}{\mathrm{e}}^{\mathrm{\Psi }}\frac{\mathit{\gamma }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }{\mathrm{)}}^{\mathrm{2}}}{\mathit{\gamma }\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}}]\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{\Psi }}\mathit{\gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}\mathit{,}\mathrm{\Psi }\right)\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }\mathit{r̂}\mathrm{2}}|}_{\mathrm{0}}\}\mathit{r̂}\mathrm{2}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}\mathit{r̂}\mathrm{4}\mathrm{\right)}\mathit{.}\end{array}$(B.12)Therefore, the models in this family are characterized by pressure isotropy in the central region, radially-biased pressure anisotropy in the intermediate part, and pressure isotropy at the boundary (since, for both and , the term of lower order in ψ is given by 2/5   ψ^5/2), at variance with the family defined by ${\mathit{f}}_{\mathrm{WT}}^{\mathrm{d}}\mathrm{\left(}\mathit{I}\mathrm{\right)}$ in which tangentially-biased anisotropy is present.

Using the asymptotic expression of the density in the central regions recorded above, in this case, the condition for the existence of the central toroidal structure is given by $\mathit{\chi }\mathit{>}\frac{\mathit{\gamma }\mathrm{(}\mathrm{1}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}}{\mathrm{12}\mathit{\gamma }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}}\left[\mathrm{3}\mathrm{+}\frac{{\mathit{C}}_{\mathrm{2}}}{\mathrm{2}}{\left(\frac{\mathrm{5}}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\right]\mathit{,}$(B.13)where C₂ is defined as in Eq. (38). By evaluating the sign of the velocity and pressure term in the radial component of the Jeans equation (see Eq. (52)), we found that, in this case, the requirement of positivity of the velocity term is just a necessary but not sufficient condition for the existence of the central toroidal structure. In other words, in this family, configurations with angular velocity higher than the angular velocity associated with the circular orbit of a single star but without a central toroidal structure can exist. The condition for maximally rotating configurations, which, in this case, is given by $\mathit{\chi }\mathit{<}\frac{\mathit{\gamma }\mathrm{(}\mathrm{1}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}}{\mathrm{8}\hspace{0.17em}\mathit{\gamma }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}\mathit{,}\mathrm{\Psi }\mathrm{)}}{\left(\frac{\mathrm{5}}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\left|}{\mathit{C}}_{\mathrm{2}}\mathrm{\right|}\mathit{,}$(B.14)completes the summary of the intrinsic properties of the family of plain-truncated differentially rotating models.

Appendix C: The numerical procedure for the construction of differentially rotating models

The construction of the equilibrium configurations for the families of differentially rotating models, defined by Eqs. (21) and (22), has been performed by numerically solving the relevant Poisson equation as a nonlinear equation for the unknown potential ψ. The code which implements the iteration procedure described in Sect. 3.2 starts with the calculation of the “seed solution”, given by a selected spherical configuration from the family of models that represent the limit in the case of vanishing internal rotation of the family of interest (i.e., defined by Eqs. (A.3) and (A.1), respectively). Such spherical solution is used, in the first step of the iteration, to evaluate (by means of a Double Gaussian Quadrature) the density distribution, given by or , on a spherical grid in the meridional plane, defined by (r_i,θ_j). The grid is linear in the radial coordinate and the angular positions are defined by ${\mathit{\theta }}_{\mathit{j}}\mathrm{=}\frac{\mathit{\pi }\mathrm{(}\mathrm{2}\mathit{j}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\hspace{0.17em}{\mathit{l}}_{\max }}\mathit{,}$(C.1)where j = 1,...,2l_max + 1. Typically, we used  ≈ 300 radial steps and we set l_max = 21, in order to have sufficient accuracy to describe the complex morphologies of the configurations in the rapid and extreme rotation regimes. The discrete direct and inverse Legendre transforms, required at each step of the iteration for the calculation of the density and potential coefficients, defined by Eqs. (26)–(25), are performed (up to the order l_max) by means of a package based on S2kit 1.0 by Kostelec & Rockmore (2004), which makes use of FFTW 3.2.1 by Frigo & Johnston (2005). The Cauchy problems for the potential radial functions expressed in Eqs. (31)–(32) are evaluated by means of a numerical integration with Romberg’s rule. The convergence condition for the solution at the n-th step of the iteration is formally defined as $\frac{{\mathit{\psi }}_{\mathit{ij}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{-}{\mathit{\psi }}_{\mathit{ij}}^{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}}{{\mathit{\psi }}_{\mathit{ij}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}}\mathit{<}{\mathrm{10}}^{-3}\mathit{,}$(C.2)for every i,j, where ${\mathit{\psi }}_{\mathit{ij}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{=}{\mathit{\psi }}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,}{\mathit{\theta }}_{\mathit{j}}\mathrm{\right)}$; about 10 (25) iteration steps are needed for the construction of configurations characterized by low (high) values of χ. The accuracy of the solutions found with our code has been checked by the following tests: (i) the virial theorem is satisfied with accuracy of the order of 10^-4 or better; (2) the radial component of the Jeans equation is satisfied with the accuracy of the order of 10^-3 or better; (3) the asymptotic behaviors, both in the central and in the outer parts, of all the moments of the distribution function are confirmed.

All Tables

All Figures