Properties of quasirelaxed stellar systems in an external tidal field
Abstract
In a previous paper, we have constructed a family of selfconsistent triaxial models of quasirelaxed stellar systems, shaped by the tidal field of the hosting galaxy, as an extension of the wellknown spherical King models. For a given tidal field, the models are characterized by two physical scales (such as total mass and central velocity dispersion) and two dimensionless parameters (the concentration parameter and the tidal strength). The most significant departure from spherical symmetry occurs when the truncation radius of the corresponding spherical King model is of the order of the tidal radius, which, for a given tidal strength, is set by the maximum concentration value admitted. For such maximally extended (or “critical”) models the outer boundary has a generally triaxial shape, given by the zerovelocity surface of the relevant Jacobi integral, which is basically independent of the concentration parameter. In turn, the external tidal field can give rise to significant global departures from spherical symmetry (as measured, for example, by the quadrupole of the mass distribution of the stellar system) only for lowconcentration models, for which the allowed maximal value of the tidal strength can be relatively high. In this paper we describe in systematic detail the intrinsic and the projected structure and kinematics of the models, covering the entire parameter space, from the case of subcritical (characterized by “underfilling” of the relevant Roche volume) to that of critical models. The intrinsic properties can be a useful starting point for numerical simulations and other investigations that require initialization of a stellar system in dynamical equilibrium. The projected properties are a key step in the direction of a comparison with observed globular clusters and other candidate stellar systems.
Subject headings:
globular clusters: general — stellar dynamics1. Introduction
It is commonly thought that globular clusters can be described as stellar systems of finite size, with a truncation in their density distribution determined by the tidal field of the hosting galaxy. Given the fact that many globular clusters undergo significant relaxation, the simplest analytical models (King 1966) that provide a successful description of these stellar systems are defined by a quasiMaxwellian distribution function , where is the singlestar energy, with a truncation prescription continuous in phase space.
Most of the interesting physical mechanisms that underlie the dynamical evolution of these stellar systems (such as evaporation and core collapse; e.g., see Spitzer 1987; Heggie & Hut 2003) depend on such truncation and are frequently studied in the context of spherical models for which the action of tides is implemented by means of the existence of a suitable truncation radius (which, in principle, for a globular cluster of given mass can be determined a priori if we know the cluster orbit and the mass distribution of the hosting galaxy) supplemented by a recipe for the escape of stars. Therefore, evolutionary models that rely on the assumption of spherical symmetry, such as Monte Carlo models (e.g., for a description of two of the codes currently used, see Giersz 1998; Joshi et al. 2000) and FokkerPlanck models (e.g., for an example of, respectively, the isotropic and anisotropic case, see Chernoff & Weinberg 1990; Takahashi et al. 1997) are necessarily based on an approximate treatment of the tidal field. In this context one interesting exception to the adoption of spherical symmetry is that of the axisymmetric FokkerPlanck models by Einsel & Spurzem (1999), which admit flattening induced by internal rotation.
Yet, if tides are indeed responsible for the truncation, they should also induce significant deviations from spherical symmetry: in the simplest case of a cluster in circular orbit about the center of the hosting galaxy, the associated (steady) tidal field is nonspherical and determines an elongation of the mass distribution in the direction of the center of mass of the hosting galaxy accompanied by a compression in the direction perpendicular to the orbit plane (e.g., see Spitzer 1987; Heggie & Hut 2003). Direct Nbody simulations, in which an external tidal field can be taken into account explicitly, provide a unique tool for the study of the evolution of a tidally perturbed cluster, especially when noncircular orbits are considered so that tidal effects are timedependent (e.g., see Baumgardt & Makino 2003); in particular, this approach has led to detailed models of the rich morphology of tidal tails, i.e. streams of stars escaped from the cluster (e.g., see Lee et al. 2006).
Indeed, deviations from spherical symmetry are observed in globular clusters (e.g., see Geyer et al. 1983; White & Shawl 1987), but, they are often ascribed to other physical ingredients, such as internal rotation. In fact, it is recognized that the issue of what determines the observed shapes of globular clusters remains unclear (e.g., see King 1961; Fall & Frenk 1985; Han & Ryden 1994; van den Bergh 2008; and references therein).
One might argue that real globular clusters are likely to be not fully relaxed, may possess some rotation and experience timedependent tides so that analytical refinements beyond the spherical onecomponent King models would not compete with the currently available numerical simulation tools (see § 5.1 for additional comments) that allow us to include these and a great variety of other detailed effects that are relevant for the quasiequilibrium configurations. However, physically simple analytical models, accompanied by the study of more realistic numerical simulations, serve as a useful tool to interpret real data and to provide insights into dynamical mechanisms, even though we know that real objects certainly include features that go well beyond such simple physical models. In this specific case, we argue that the triaxial geometry is a natural attribute of the physical picture of tidal truncation, which has already proved to provide useful guidance in the study of globular clusters.
As demonstrated in a previous paper (Bertin & Varri 2008; hereafter Paper I), a complete characterization of the physically simple description of globular clusters as quasirelaxed tidallytruncated stellar systems, in terms of fully selfconsistent triaxial models, can be provided. In order to better address the role of tides in determining the observed structure of globular clusters, in Paper I we have thus constructed selfconsistent nonspherical equilibrium models of quasirelaxed stellar systems, obtained from the spherical case by including in their distribution function the effects associated with the presence of an external tidal field explicitly. We recall that our models consider the stellar system in circular orbit within the hosting galaxy, for simplicity assumed to be sphericallysymmetric. Therefore, in the corotating frame of reference, the tidal field experienced by the system is static and the Jacobi integral is available. In this physical picture the typical dynamical time associated with the orbits inside the cluster is assumed to be much smaller than the external orbital time. The procedure described in the previous paper starts by replacing the singlestar energy with the Jacobi integral in the relevant distribution function if , with the cutoff constant, and otherwise. Thus the collisionless Boltzmann equation is satisfied. The construction of the selfconsistent models then requires the solution of the associated PoissonLaplace equation, that is of a secondorder elliptic partial differential equation in a free boundary problem, because the boundary of the configuration, which represents the separation between the Poisson and Laplace domains, can be determined only a posteriori. The idea of using the Jacobi integral for the construction of tidal triaxial models had been proposed also by Weinberg (1993). A family of models based on the distribution function , but constructed with a different method, was discussed by Heggie & Ramamani (1995); a comparison with our models will be given in § 3.2.
Here we describe the properties of the resulting twoparameter family of physically justified triaxial models, constructed with the method of matched asymptotic expansions (see § 4 and § A in Paper I), inspired by previous studies of the formally similar problem of rotating polytropes (Smith 1975; 1976). In particular, we illustrate the properties of the first and secondorder solutions.
The paper is organized as follows. In § 2 we present a thorough description of the relevant parameter space. The intrinsic and projected density distributions are discussed in § 3, with special emphasis on the global and local quantities that can be used as diagnostics of deviations from spherical symmetry. Intrinsic and projected kinematics are addressed in § 4. The concluding § 5 gives the summary of the paper with a discussion of the results obtained and a comment on the complex physical phenomena that a large body of evolutionary models based on numerical investigations has shown to characterize the periphery of globular clusters.
2. The Parameter Space
The triaxial tidal models are characterized by two physical scales (corresponding to the two free constants and in the distribution function ) and two dimensionless parameters. The latter parameters are best introduced by referring to the formulation of the Poisson equation in terms of the dimensionless escape energy
(1) 
where is the dimensionless cluster meanfield potential (to be determined selfconsistently) and represents the tidal potential (with the numerical coefficient , where and are respectively the epicyclic and orbital frequency, depending on the potential of the hosting galaxy); the hat on the spatial coordinates denotes that they are measured in units of the scale length .
Then the first parameter, already available in the spherical case, is the concentration of the system and can be expressed as a dimensionless measure of the central depth of the potential well: . The second parameter, the tidal strength in Eq. (1), is defined as
(2) 
i.e. as the ratio of the square of the orbital frequency (of the revolution of the stellar system around the center of the hosting galaxy) to the square of the dynamical frequency associated with the central density of the stellar system. Alternatively, the effect of the tidal field can be measured by the extension parameter
(3) 
where is the truncation radius of the spherical King model characterized by the same value of and is the tidal (or Jacobi) radius, i.e. the distance from the origin (the center of the stellar system) of the two nearby Lagrangian points of the restricted threebody problem considered in our simple physical picture. A given model will be labelled by the pair of values or, equivalently, by the pair . The dimensionless cutoff constant can be expressed in a natural way as an asymptotic series with respect to the tidal parameter where the terms , as discussed in Paper I, depend only on .
Much like the Hill surfaces for the standard restricted threebody problem, we now consider the family of zerovelocity surfaces defined by the condition , which represents the boundary of our models. These surfaces can be open or closed, depending on the value of the cutoff constant , which is determined by the selected values of the two dimensionless parameters that characterize the model. To be consistent with the hypothesis of stationarity, we only consider closed configurations. We call “critical models” those that are bounded by the critical zerovelocity surface (which is the outermost available closed surface). For each value of , the critical value of the tidal parameter can be found by (numerically) solving the system
(4) 
where the unknowns are and . The method of matched asymptotic expansions proposed in Paper I for the solution of the relevant PoissonLaplace equation requires an expansion in spherical harmonics, therefore it can be easily recognized that the first condition of Eq. (4) is equivalent to the requirement of vanishing gradient, which identifies the saddle points of the critical surface. In the general case, the condition determines the value of for a given tidal strength , therefore .
By using in the escape energy defined in Eq. (1)
the zerothorder expression for the cluster potential and for the cutoff constant
(5) 
thus leading to an expression for in terms of the dimensionless truncation radius of the corresponding spherical King model and a first estimate of the critical value of
(6) 
(7) 
The first expression can also be written as (see Spitzer 1987). We recall that and that . Therefore, the righthand side of Eqs. (6) and (7) depends only on the value of .
Here one important comment is in order. Strictly speaking, the complete solution for derived by the method of matched asymptotic expansions in § 4 of Paper I is a welljustified global uniform solution only for subcritical (underfilled) models. Close to the condition of criticality, i.e., when , in the vicinity of the boundary surface the tidal term (which is considered a small correction in the construction of the asymptotic solution) becomes comparable to the cluster term , so that the asymptotic solution is expected to break down. For such models, the iteration method described in § 5.2 of Paper I, which does not rely on the assumption that the tidal term is small, is preferred and expected to lead quickly to more accurate solutions. In practice, in line with previous work on the similar problem for rotating polytropes mentioned in the Introduction, we argue that the use of the secondorder asymptotic solution constructed in Paper I will give sufficiently accurate solutions in the determination of the critical value of the tidal parameter from the system in Eq. (4) and in the consequent assessment of the general properties of models even when close to the critical case. The main reason at the basis of this argument is that, even for closetocritical models, only very few stars populate the region where the asymptotic analysis breaks down, so that the overall solution should be only little affected. A direct comparison between selected critical models calculated with both the perturbation and the iteration method is presented in Appendix A.
If we make use of the full secondorder asymptotic solution for the escape energy , in which the secondorder expressions for the cluster potential (recorded in Eq. (B3)) and the cutoff constant are used, the system in Eq. (4) can be rearranged and written in standard form
(8) 
with
(9) 
(10) 
(11) 
(12) 
(13) 
(14) 
where the relevant constants ( with , with and with ), which are determined by the matching process, are defined in § 4.4 of Paper I. [The corresponding system based on the firstorder solution for can be recovered by setting .] This system has been solved numerically, by means of the NewtonRaphson method, since the equations are nonlinear in (in particular, they are polynomials of fifth and seventhorder, for the first and secondorder solution respectively). As noted in the discussion of the simpler Eq. (5), the solution for can then be represented as a function of the concentration parameter . The parameter space of the first and secondorder models has been explored by means of an equally spaced grid from to at steps of and from to at steps of .
The parameter space for the secondorder models is presented in Fig. 1. The plot provides the contour levels of the extension parameter , with the uppermost solid line corresponding to (thus identifying the critical models), based on the choice (Keplerian host galaxy). The critical curves for (host galaxy characterized by flat rotation curve) and for (Plummer potential evaluated at , with the model scale radius) are shown as a dashed and dotdashed line, respectively.
Subcritical, underfilled models (bottomleft corner of the figure), with , are only little affected by the tidal perturbation. The maximally deformed models are those with (close to the uppermost solid line, i.e. closetocritical configurations). Figure 1 shows that the critical value for the tidal strength parameter depends strongly on concentration, with a variation of almost four orders of magnitude in the explored range of . The figure also indicates that for lower values of the critical curve moves upwards, i.e. the available parameter space increases.
The difference between the critical value of the tidal strength parameter for first and secondorder models (for a chosen value of ) is very small, around for lowconcentration models, down to or less for models with . The critical value of the extension parameter depends only weakly on concentration and on , as illustrated in Fig. 2.
In closing this section, se should reiterate that, in spite of the abundant use of symbols required by the analysis, the family of models that we have studied is characterized by two dimensionless parameters . [As an alternative pair, we may refer to the standard concentration parameter , equivalent to and frequently used in the context of spherical King models, and to the extension parameter , equivalent to .] The free constants and that appear in the distribution function set the two physical scales. In turn, the models constructed by Heggie & Ramamani (1995) are a oneparameter family of models, because these authors focused on the critical case and did not discuss the subcritical regime.
3. Intrinsic and projected density distribution
3.1. Intrinsic density profile
The models are characterized by reflection symmetry with respect to the three natural coordinate planes. With respect to the unperturbed configuration (i.e. the spherical King model with the same value of ), they exhibit an elongation along the axis (defined by the direction of the center of the host galaxy), a compression along the axis (the direction perpendicular to the orbit plane of the globular cluster), and only a very modest compression along the axis.
Models with , regardless of the value of , are practically indistinguishable from the corresponding spherical King models; significant departures from spherical symmetry occur for models with or higher. In Fig. 3.a we show the density profile along the axis for a selection of critical secondorder models with in comparison with that of the corresponding spherical King models; note that for a model with the elongation is already significant at , while for a model with a similar elongation is reached only at much lower density levels (). The corresponding profiles along the axis and the axis are given in Fig. 3.b.
For completeness, we checked the dependence of our density profiles on the potential of the host galaxy. Consistent with the general trends suggested by Fig. 2, the elongation along the axis and the compression along the axis for the models with turn out to be slightly weaker than for the models with smaller values of .
Since the dimensionless density distribution of a model, identified by , is given by , where is the dimensionless escape energy and is a monotonically increasing function that vanishes for vanishing argument (see Eq. (11) in Paper I), there is a onetoone correspondence between isodensity and isovelocity surfaces, the latter being defined by the condition , where is a constant (with ). We recall that nonspherical models often exhibit equipotential surfaces rounder than the isodensity surfaces (e.g., see Evans 1993; Ciotti & Bertin 2005). The reason for the presence of this property in our models is that the supporting distribution function depends only on the Jacobi integral, i.e. the (isolating) energy integral in the rotating frame. Therefore, for each value of , we can define the semiaxes , , and of the corresponding triaxial isodensity surface, by means of the intersections of the surface with the , , and axes, which turn out to follow the ordering . The shape of the triaxial configuration can thus be described in terms of the polar and equatorial eccentricities, defined as and , respectively.
A surprising result can be derived analytically. In the innermost region (i.e., for ), the dimensionless escape energy can be expanded to second order in the dimensionless radius
(15)  
Here some terms of the secondorder solution do not contribute (e.g., it can be readily checked that and ). Then by setting , we find that in the innermost region the eccentricities tend to the following nonvanishing central values
(16) 
(17) 
where
(18) 
(19) 
which depend explicitly on the tidal strength and implicitly on the concentration. This result is nontrivial. In fact, since the tidal potential is a homogeneous function of the spatial coordinates, naively we might expect that in their central region the models reduce to a perfectly spherical shape (i.e., ), even for finite values of the tidal strength. Instead, and are and strictly vanish only in the limit of vanishing tidal strength.
Figure 4 shows the central values of the eccentricities for secondorder models with and selected values of concentration, as a function of tidal strength within the range . Consistent with the general trends identified in the discussion of the parameter space, lowconcentration models show the most significant departures from spherical symmetry. The full eccentricity profiles (as a function of the major axis) are shown in Fig. 5 for a selection of critical secondorder models; here the calculation of and has been performed by numerically determining the values of the semiaxes of a number of reference isovelocity surfaces, defined by with . Outside the central region, the profiles increase monotonically and, independently of concentration, at the boundary they reach approximately a fixed maximum value ( and ), which corresponds to the fact that the shape of the boundary surface of a critical model (; see also Fig. 2) depends only modestly on concentration.
3.2. Comparison with the models constructed by Heggie & Ramamani (1995)
The method used in Paper I for the construction of the models illustrated in this paper can be summarized as follows. The solution in the internal (Poisson) and external (Laplace) domains are expressed as an asymptotic series with respect to the dimensionless parameter , representing the tidal strength (defined in Eq. (2)), which is considered to be small. The thorder term of the asymptotic series of the internal and external solution are denoted by and respectively, so that the zerothorder terms define the standard spherical King models (see § 4.1 in Paper I for a detailed discussion). The quantity indicates the thorder external solution, i.e. the corresponding asymptotic series truncated at the term ; a similar notation holds for the internal solution. The validity of the expansion breaks down where the second term is comparable to the first, i.e. where . This singularity is cured by introducing a boundary layer in which both the spatial coordinates and the solution are suitably rescaled with respect to the tidal parameter. To obtain a uniformly valid solution over the entire space, an asymptotic matching (see Van Dyke 1975; Eq. (5.24)) is performed between the pairs and . Each term is then expanded in spherical harmonics with radial coefficients . The internal region requires a numerical solution of the Cauchy problems for the radial coefficients (we used a fourthorder RungeKutta code) while in the external region a formal solution with multipolar structure is available and in the boundary layer the integration in the radial variable can be performed analytically (see § 4.2 and 4.3 in Paper I, respectively).
The models described by Heggie & Ramamani (1995) are also based on a perturbation approach, but the method used is different from ours, provides a solution of the Poisson equation that is firstorder with respect to the tidal parameter, and is restricted to the “critical” case; actually, as noted in § 2 after Eq. (7), the perturbation approach is bound to break down in the critical case. Technically speaking, their method is in the form of a “patching” procedure, in contrast with our asymptotic matching. Therefore, the models constructed in Paper I, while consistent, to first order, with those of Heggie & Ramamani (1995), are more general. We also recall that our method is also applicable to systems described by different distribution functions (see Appendix B and C in Paper I).
We have thus performed a quantitative comparison between the
intrinsic density profiles of our critical secondorder models and those of
the models
3.3. Global quantities
The previous discussion has focused on the shape of the isodensity surfaces of the models. In particular, some interesting conclusions have been derived based on a local analysis of the central region and of the outer boundary of the configuration. The maximal departures from spherical symmetry are reached at the periphery, but these hinge on the distribution of the very small number of stars that populate the outer region of the cluster. We may thus wish to study some global quantities that better characterize whether significant amounts of mass (and, correspondingly, of light) are actually involved in the deviation of the model from spherical symmetry. One standard such global measure is provided by the quadrupole moment tensor
(20) 
with the integration to be performed in the volume of the entire configuration. Here the notation for the function and for the constant is the same as in Eq. (11) of Paper I. In contrast with the frequently used inertia tensor (e.g., see Chandrasekhar 1969; chap. 2), the quadrupole moment is defined in such a way that in the spherical limit it vanishes identically. In our coordinate system it is diagonal. Note that the tidal distortions require that the nonvanishing terms of the inertia tensor follow the ordering ; to visualize the geometry of the system, we may thus also refer to the average polar and equatorial eccentricities and defined by the relations , . In general, we have and , with the prolate configuration identified by , i.e. .
Since most of the mass is contained in the inner regions, global quantities can be evaluated approximately by neglecting the contribution from the region corresponding to the boundary layer. We can thus use the secondorder solutions for obtained in Paper I by the method of matched asymptotic expansions and conveniently reduce the calculation of global quantities to an easier integration in spherical coordinates inside the sphere of radius . Therefore, for the quadrupole moment tensor we find
(21) 
We emphasize that this estimate is expected to be a good approximation only for those models for which the contribution of the boundary layer is negligible with respect to the one of the internal sphere of radius .
The relevant components on the diagonal can be expressed in terms of the matching constants of the external solution (for the relevant definitions, see Eqs. (63) and (69) in Paper I)
(22) 
(23) 
(24) 
We recall that the constants and are positive, while and are negative (and larger in magnitude). Therefore, is positive and and are negative, consistent with the detailed elongation and compressions observed in the density profile. A summary of the derivation of these formulae is provided in Appendix B.
As a measure of the degree of triaxiality of a given configuration, we have calculated the following ratios
(25) 
(26) 
which depend explicitly on the tidal strength parameter and implicitly on the concentration parameter. In the limit of vanishing tidal strength, we find
(27) 
(28) 
where are the quadrupole coefficients of the tidal potential (see Eqs. (38) and (39) in Paper I). This result is nontrivial, because, in this limit, numerator and denominator are both expected to vanish. Note that only for the ratio .
Earlier in this paper we mentioned that two physical scales, such as total mass and central velocity dispersion, correspond to the two dimensional constants and that appear in the distribution function . In fact, the total mass of the system is given by
(29) 
If we insert the secondorder solution for obtained in Paper I, we find
(30) 
with
(31) 
Here each term of the expansion is related to the corresponding constant with (i.e., the monopole term) of the expression of the external solution of the PoissonLaplace equation calculated by means of the method of matched asymptotic expansion (for the relevant definitions, see Eqs. (59),(61),(67) in Paper I).
The quality of the analytical estimates for the total mass and the quadrupole moment tensor has been checked by comparing the values obtained from asymptotic analysis with those resulting from direct numerical integration of Eq. (29) and Eq. (3.3) respectively, in which the density profile is used without any additional expansion. The integration of the distribution function over the entire space, required by those global quantities, has been performed by means of an Adaptive Monte Carlo method (the algorithm VEGAS, see Press et al. 1992; § 7.8), well suited for our geometry. For the quadrupole, the results are illustrated in Fig. 7. For the mass, the dimensionless function is basically unchanged (within 0.5 percent) with respect to the function characterizing the spherical King models; the Monte Carlo integration is very accurate, with relative errors around , and the analytical approximation given by Eq. (3.3) shows an excellent agreement for every value of in the whole range of the extension parameter .
The average eccentricities for critical secondorder models as a function of concentration, with , for two different choices of the host galaxy potential (), are shown in Fig. 8. The values are calculated directly from the definitions given earlier in this subsection in terms of with no additional expansions. A nonmonotonic dependence on concentration is revealed, with generally higher average eccentricities attained by lowconcentration models. The trends of the polar and equatorial eccentricities are basically the same, as shown by the fact that the related curves in the plot can be matched approximately by a rigid translation. As expected (see § 3.1), models with show a larger separation between polar and equatorial eccentricities than models with . The presence of a minimum for the curves at occurs approximately at the location where the function shows an inflection point (regardless of the value of ; see Fig. 1).
3.4. Projected density profile
Under the assumption of a constant masstolight ratio, projected models can be compared with the observations. We have then computed surface (projected) density profiles and projected isophotes. The projection has been performed along selected directions, identified by the viewing angle , corresponding to the axis of a new coordinate system related to the intrinsic system by the transformation ; the rotation matrix is expressed in terms of the viewing angles, by taking the axis as the line of nodes (see Ryden 1991; for an equivalent projection rule). Given the symmetry of our models (see § 3.1), viewing angles can be chosen from the first octant only. In particular, we used a polar grid defined by and with , and calculated (by numerical integration, using the Simpson rule) the dimensionless projected density
(32) 
where with the edge of the cluster along the axis of the intrinsic coordinate system (i.e., we “embedded” the triaxial configuration in a sphere of radius given by its maximal extension). The projection plane has been sampled on an equallyspaced square grid centered at the origin (note that for the projected density is correctly set to zero).
The morphology of the isophotes of a given projected image can be described in terms of the ellipticity profile, defined as where and are the semiaxes, as a function of the major axis . As already noted for the (intrinsic) eccentricity profiles, the deviation from circularity increases with the distance from the origin. In the inner region, the ellipticity is consistent with the central eccentricities and calculated in § 3.2.
By taking lines of sight different from the axes of the symmetry planes, we have also checked whether the projected models would exhibit isophotal twisting. For all the cases considered, the position angle of the major axis remains unchanged over the entire projected image. Tests made by changing the resolution of the grid confirm that, even in the most triaxial case (), no significant twisting is present.
The first two panels of Fig. 9 show the projected images of a critical secondorder model with and along the and directions, (i.e., the and axis of the intrinsic system), corresponding, respectively, to the least and to the most favorable line of sight for the detection of the intrinsic flattening of the model. For the same model, the third panel illustrates the ellipticity profiles for various lines of sight.
Figure 10 shows the surface density profiles along the and axes of the projection plane for ten critical secondorder models with , viewed along the direction. As a further characterization, for the same models in the lower panel we also present the surface density profiles obtained by averaging the projected density distribution on circular annuli; this conforms to the procedure often adopted by observers in dealing with density distributions with very small departures from circular symmetry (e.g., see Lanzoni et al. 2007). As expected, circularaveraged profiles lie between the corresponding regular profiles taken along the principal axes of the projected image.
4. Intrinsic and projected kinematics
By construction, the models are characterized by isotropic velocity dispersion. The intrinsic velocity dispersion can be determined directly as the second moment in velocity space (normalized to the intrinsic density) of the distribution function
(33) 
where represents the incomplete gamma function (near the boundary of the configuration, the velocity dispersion profile scales as ). This shows that the isodensity surfaces of the models are in a onetoone correspondence with the isovelocity and isobaric surfaces (defined by const). As noted for the intrinsic density profiles in § 3.1, a compression along the axis and an elongation along axis occur also for the intrinsic velocity dispersion profiles. In Fig. 11.a we present the intrinsic velocity dispersion profiles along the axis for the same critical models illustrated in Fig. 3 compared to the profiles of the corresponding spherical King models. The behavior of the projected velocity dispersion profiles near the boundary is significantly different from that of the spherical models.
The projected velocity moments can be calculated by integrating along the line of sight (weighted by the intrinsic density) the corresponding intrinsic quantities. Therefore, the projected velocity dispersion is given by
(34) 
Figure 11.b shows the projected velocity dispersion profiles along the and axis of the projection plane for the same models displayed in Fig. 10 (the line of sight is defined by ).
5. Discussions and conclusions
5.1. A comment on the complex structure of the outer regions
In Paper I we have noted the singular character of the mathematical problem associated with the free boundary set by the threedimensional surface of these tidally truncated models. In § 2 of the present Paper, we have further emphasized the additional singularity that characterizes closetocritical models (see comment after Eq. (7), which prompted the analysis described in the Appendix A). As is true in general in the study of boundary layers and similar problems, it is no wonder that in the vicinity of such critical boundaries, a number of complex physical effects may take place and play an important role in determining the detailed structure of the solution. On the other hand, the properties of the derived solution away from the boundary are quite robust (see § 3.3 and the Appendix). As to some of the subtle properties of the expected distribution function close to the edge of a cluster, it is interesting to summarize here below the main results that emerge from a vast body of evolutionary models (Nbody, FokkerPlanck, Monte Carlo, gas), on the issue of the interplay between pressure anisotropy and tidal effects.
Since the first solutions of the FokkerPlanck equation by means of a Monte Carlo approach, as described by Hénon (1971) or Spitzer & Hart (1971), it has been shown that onecomponent isolated spherical clusters, starting from a variety of initial conditions (see Spitzer 1987; chap.4 and references therein), develop a corehalo structure in which the central parts are almost isotropic while the outer regions are characterized by radial anisotropy. A commonly reported interpretation is that the halo is mainly populated by stars scattered out from the core on radial orbits. If the evolution of a cluster is influenced not only by internal processes but also by the presence of an external tidal field, the growth of pressure anisotropy is significantly damped. Direct Nbody simulations (e.g., see Giersz & Heggie 1997; Aarseth & Heggie 1998; Baumgardt & Makino 2003; Lee et al. 2006), anisotropic FokkerPlanck (e.g., see Takahashi et al. 1997), and Monte Carlo models (e.g., see Giersz 2001) of both single and multimass systems suggest that clusters in circular orbits (and even in eccentric orbits, see Baumgardt & Makino 2003), starting from isotropic initial conditions, tend to preserve pressure isotropy, except for the outermost parts which become tangentially anisotropic due to the preferential loss of stars on radial orbits, induced by the tidal field. The overall agreement on this result is nontrivial, because of the aforementioned differences in the treatment of the external tidal field. Even extreme cases of timedependent tides, such as disk shocking, influence the degree of pressure anisotropy since it has been shown that they may represent a dominant mechanism (“shock relaxation”) of the energy redistribution, leading to a substantial isotropy, of the weakly bound stars (see Oh & Lin 1992; Kundic & Ostriker 1995; both papers are based on a FokkerPlanck approach).
These results confirm that, of course, the properties of the models constructed in Paper I and described in this paper should be taken only as a zerothorder reference frame, to single out the precise role of external tides, and should not be taken literally as a realistic representation of real objects since a number of simplifying assumptions are made. On the other hand, by comparing data with such an idealized reference model it will be possible to better assess the role of tides with respect to other physical ingredients studied separately.
5.2. Summary and concluding remarks
The main results that we have obtained from the detailed analysis of the family of tidal triaxial models can be summarized as follows:

Two tidal regimes exist, namely of low and highdeformation, which are determined by the combined effect of the tidal strength of the field and of the concentration of the cluster. The degree of deformation increases with the degree of filling of the relevant Roche volume. Far from the condition of Roche volume filling, the models are almost indistinguishable from the corresponding spherical King models. A number of studies have investigated the evolution of tidally perturbed stellar systems initially underfilling their Roche lobe (e.g., see Gieles & Baumgardt 2008; Vesperini et al. 2009), concluding that some of the relevant dynamical processes, in particular evaporation, depend on the degree of filling of the Roche volume. The intrinsic properties of the models discussed in this paper can be useful for setting selfconsistent nonspherical initial conditions of numerical simulations aiming at studying in further details the evolution of configurations starting from subcritical tidal equilibria.

For a given tidal strength, there exists a maximum value of concentration for which a closed configuration is allowed (see Fig. 1). In such “critical” case, the truncation radius of the corresponding spherical King models is of the same order of the tidal radius of the triaxial model. The shape of the boundary of the maximally deformed models is given by the “critical” zerovelocity surface of the relevant Jacobi integral and is basically independent of concentration, while the deformation of the internal region strongly depends on the value of the critical tidal strength and is more significant for lowconcentration models. This statement agrees with a general trend noted by White & Shawl (1987) for the globular cluster system of our Galaxy.

The structure of the models can be described in terms of the polar and equatorial eccentricity profiles of the intrinsic isodensity surfaces. The maximal departures from spherical symmetry are reached in the outskirts.

For finite tidal strength, the central values of the polar and equatorial eccentricities are finite, ; this result is nontrivial since the tidal potential which induces the perturbation is a quadratic homogeneous function of the spatial coordinates.

Global measures of the degree of triaxiality in terms of the quadrupole moment tensor have been introduced and calculated for different values of the tidal strength and different potentials of the host galaxy. The potential of the host galaxy sets the geometry of the tidal perturbation, as nicely shown by the analytic expressions for the ratios of the components of the quadrupole moment tensor, which, in the case of firstorder models, reduce to simple functions of the coefficient . Average eccentricities have been calculated from the inertia tensor components, evaluated numerically.

As a key step in the direction of a comparison with observations, projected density profiles and ellipticity profiles have been calculated for a number of models for several lines of sight.

The study of the relevant (projected) isophotes indicates that no isophotal twisting occurs. This result is nontrivial since the models are nonstratified and centrallyconcentrated, conditions under which, in principle, isophotal twist may occur (see Stark (1977); models based on Stäckel potentials are also known to be twistfree, as shown by Franx (1988)).

Finally, close to the boundary, the intrinsic and projected kinematics shows significant differences with respect to that of spherical models.
Since our models are all characterized by monotonically increasing ellipticity profiles, they cannot explain the variety of behavior of observed ellipticity profiles (see Geyer et al. 1983), but the range of the predicted flattening () is consistent with that observed in most globular clusters (see White & Shawl (1987) for the clusters in the Milky Way and Frenk & Fall (1982) for those in the LMC). Therefore, with the exception of special clusters such as Omega Centauri, we think that the modest but frequently observed deviations from spherical symmetry might have their origin traced to tides.
Finally, since our models are intrinsically more elongated (in the direction of the center of the host galaxy) than spherical King models, they might be useful for interpreting clusters with a surface brightness profile extending beyond that predicted by the spherical King models. Recent investigations (see McLaughlin & van der Marel 2005) suggest that such “extratidal” structures are a fairly generic feature, especially for extragalactic clusters, and not just a transient property, present only at young ages.
On the basis of the work presented here, in a following paper we plan to address in detail the issues raised by observations.
Appendix A Perturbation vs. iteration
For completeness, we calculated selected models also by means of the iteration method described in § 5.2 of Paper I, in order to verify the quality of the solution obtained with the method of matched asymptotic expansions, in particular in the critical regime. This technique follows the approach proposed by Prendergast & Tomer (1970). The basic idea is to get an improved solution of the Poisson equation (see Eq. (77) in Paper I) by evaluating the righthand side with the solution obtained from the immediately previous step, starting from the spherical King models taken as initial “seed solutions”. In our code, the required spherical harmonic analysis and synthesis of density and potential have been performed by means of the S2kit 1.0 package (Kostelec & Rockmore 2004), which makes use of FFTW 3.2.1 (Frigo & Johnson 2005). We decided to truncate the harmonic series at in order to be consistent with the maximum harmonic index admitted by the secondorder asymptotic solution. The iteration stops when convergence to four significant digits in the whole domain of the solution is reached.
For selected values of the concentration parameter in the range (for simplicity, we considered only the case of an external potential with ), the corresponding critical value of the tidal strength parameter obtained with the iteration method is consistent to with the value determined by the numerical solution of Eq. (8), in which the constants obtained from the asymptotic matching are used. For a critical model up to 20 iteration steps are required for convergence, while a subcritical solution typically takes in 4 to 8 steps.
Appendix B Global quantities from the multipole expansion of the cluster potential
Based on the expansion in spherical harmonics of given in Eq. (3.70) of Jackson (1999), the external potential generated by our model can be expressed by means of the following multipole expansion
(B1) 
where the multipole coefficients are defined as
(B2) 
with the integration to be performed in the volume with boundary surface defined by . This general expression can be compared with the secondorder solution of the Laplace equation obtained in Paper I
(B3)  
in order to determine the relation between the secondorder multipole coefficients (i.e. calculated by means of the secondorder expression for the density) and the matching constants that appear in Eq. (B3). Therefore, we find that