QUASI-RELAXED STELLAR SYSTEMS IN TIDAL FIELD

# Properties of quasi-relaxed stellar systems in an external tidal field

## Abstract

In a previous paper, we have constructed a family of self-consistent triaxial models of quasi-relaxed stellar systems, shaped by the tidal field of the hosting galaxy, as an extension of the well-known 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 zero-velocity 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 low-concentration 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 sub-critical (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.

globular clusters: general — stellar dynamics

## 1. 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 quasi-Maxwellian distribution function , where is the single-star 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 Fokker-Planck 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 Fokker-Planck 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 N-body 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 non-circular orbits are considered so that tidal effects are time-dependent (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 time-dependent tides so that analytical refinements beyond the spherical one-component 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 quasi-equilibrium 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 quasi-relaxed tidally-truncated stellar systems, in terms of fully self-consistent 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 self-consistent non-spherical equilibrium models of quasi-relaxed 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 spherically-symmetric. 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 single-star energy with the Jacobi integral in the relevant distribution function if , with the cut-off constant, and otherwise. Thus the collisionless Boltzmann equation is satisfied. The construction of the self-consistent models then requires the solution of the associated Poisson-Laplace equation, that is of a second-order 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 two-parameter 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 second-order 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

 ψ(^r)=aH0−[aΦC(^r)+ϵT(^x,^z)] , (1)

where is the dimensionless cluster mean-field potential (to be determined self-consistently) 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

 ϵ≡Ω24πGρ0 , (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

 δ≡^rtr/^rT , (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 three-body 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 cut-off 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 three-body problem, we now consider the family of zero-velocity 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 cut-off 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 zero-velocity 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

 {∂^xψ(^x=^rT,^y=0,^z=0;ϵcr)=0ψ(^x=^rT,^y=0,^z=0;ϵcr)=0 , (4)

where the unknowns are and . The method of matched asymptotic expansions proposed in Paper I for the solution of the relevant Poisson-Laplace 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 zeroth-order expression for the cluster potential and for the cut-off constant1, the system in Eq. (4) becomes

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩λ0^r2T+9ϵcrν^rT=0α0−λ0^rT+92ϵcrν^r2T=0 , (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

 ^r(0)T=32λ0α0=32^rtr , (6)
 ϵ(0)cr=−8α30243νλ20=−8243ν1^r3trλ0 . (7)

The first expression can also be written as (see Spitzer 1987). We recall that and that . Therefore, the right-hand 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 well-justified global uniform solution only for sub-critical (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 second-order 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 close-to-critical 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 second-order asymptotic solution for the escape energy , in which the second-order expressions for the cluster potential (recorded in Eq. (B3)) and the cut-off constant are used, the system in Eq. (4) can be re-arranged and written in standard form

 {A(^rT)ϵ2cr+2B(^rT)ϵcr+C(^rT)=0D(^rT)ϵ2cr+2E(^rT)ϵcr+F(^rT)=0 , (8)

with

 A(^rT)=λ2^r4T−(b20−b22√3)34√5π^r2T+b404516√π −b42158√5π+b441516√35π , (9)
 B(^rT)=λ1^r4T−(a20−a22√3)34√5π^r2T+9ν^r7T , (10)
 C(^rT)=2λ0^r4T , (11)
 D(^rT)=α2^r5T−λ2^r4T+(b20−b22√3)14√5π^r2T −b40916√π+b4238√5π−b44316√35π , (12)
 E(^rT)=α1^r5T−λ1^r4T+(a20−a22√3)14√5π^r2T+92ν^r7T , (13)
 F(^rT)=2(α0^r5T−λ0^r4T) , (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 first-order solution for can be recovered by setting .] This system has been solved numerically, by means of the Newton-Raphson method, since the equations are nonlinear in (in particular, they are polynomials of fifth and seventh-order, for the first and second-order 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 second-order 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 second-order 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 dot-dashed line, respectively.

Sub-critical, underfilled models (bottom-left 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. close-to-critical 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 second-order models (for a chosen value of ) is very small, around for low-concentration 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 one-parameter family of models, because these authors focused on the critical case and did not discuss the sub-critical 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 second-order 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 one-to-one 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 semi-axes , , 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

 ψ(int)(^r)∼Ψ−32^r2+ϵ^r2[−32(1−ν)+A20Y20(θ,ϕ)+ (15) A22Y22(θ,ϕ)]+ϵ22^r2[(1+B20)Y20(θ,ϕ)+(1+B22)Y22(θ,ϕ)] .

Here some terms of the second-order 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 non-vanishing central values

 e0={ϵ(A22−√3A20+(ϵ/2)[1+B22−√3(1+B20)])√15/π}1/2{6+2ϵ[3(1−ν)−A20√5/π]−ϵ2(1+B20)√5/π}1/2 , (16)
 η0={ϵ[2A22+(1+B22)ϵ]√15/π}1/2{6+ϵd1+(ϵ2/2)d2}1/2 , (17)

where

 d1=6(1−ν)+(A20+√3A22)√5/π (18)
 d2=[1+B20+√3(1+B22)]√5/π (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 second-order 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, low-concentration 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 second-order models; here the calculation of and has been performed by numerically determining the values of the semi-axes 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 th-order term of the asymptotic series of the internal and external solution are denoted by and respectively, so that the zeroth-order terms define the standard spherical King models (see § 4.1 in Paper I for a detailed discussion). The quantity indicates the th-order 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 fourth-order Runge-Kutta 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 first-order 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 second-order models and those of the models2 by Heggie & Ramamani (1995), both referred to the case in which the host galaxy is Keplerian (). As desired, there is substantial consistency except for the outermost part of the boundary layer in which our models are slightly more compact, due to a global “boxiness” effect induced by the second-order term present in our models in which harmonics of order also play a role. In Fig. 6 we represent the relative difference between the two density profiles evaluated along the -axis (with the coordinates scaled with respect to the truncation radius instead of the usual scale radius ) for selected values of . A similar behaviour is found also along the -axis, for in the range , and along the -axis, for in the range . In the central part of the internal region, along the principal axes, the relative difference is smaller than 5 percent for every value of we tested, while near the transition to the boundary layer (i.e and ) a difference of 20 percent can be reached in the case of highly concentrated models (). We interpret these differences as due to the combined effects of the patching vs. matching adopted process and of the different grid on which the Cauchy problems for the radial coefficients are solved (we used a regular radial grid while Heggie & Ramamani (1995) used a more complex tabulation resulting from their choice of taking the zeroth-order cluster potential as the independent variable and of the function instead of ).

### 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

 Qij=∫V(3xixj−r2δij)ρ(r)d3r= ^Ar50∫V(3^xi^xj−^r2δij)^ρ(^r)d3^r=^Ar50^Qij , (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 non-vanishing 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 second-order 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

 ^Q(2)ij=^Qij,1ϵ+^Qij,2ϵ22 . (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)

 ^Q(2)xx=29√5π^ρ(Ψ)[(a20ϵ+b20ϵ22) −√3(a22ϵ+b22ϵ22)] , (22)
 ^Q(2)yy=29√5π^ρ(Ψ)[(a20ϵ+b20ϵ22) +√3(a22ϵ+b22ϵ22)] , (23)
 ^Q(2)zz=−49√5π^ρ(Ψ)(a20ϵ+b20ϵ22) . (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

 ^Q(2)yy^Q(2)xx=(a20+b20ϵ/2)+√3(a22+b22ϵ/2)(a20+b20ϵ/2)−√3(a22+b22ϵ/2) , (25)
 ^Q(2)zz^Q(2)xx=−2(a20+b20ϵ/2)(a20+b20ϵ/2)−√3(a22+b22ϵ/2) , (26)

which depend explicitly on the tidal strength parameter and implicitly on the concentration parameter. In the limit of vanishing tidal strength, we find

 ^Qyy^Qxx∼^Q(1)yy^Q(1)xx=T20(^rtr)+√3T22(^rtr)T20(^rtr)−√3T22(^rtr)=−ν−12ν+1 , (27)
 ^Qzz^Qxx∼^Q(1)zz^Q(1)xx=−2T20(^rtr)T20(^rtr)−√3T22(^rtr)=−2+ν2ν+1 , (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

 M=∫Vρ(r)d3r=^Ar30∫V^ρ(^r)d3^r=^Ar30^M . (29)

If we insert the second-order solution for obtained in Paper I, we find

 ^M(2)=∫2π0dϕ∫π0dθsinθ∫^rtr0d^r^r2^ρ(2)(^r,θ,ϕ) =^M0+ϵ^M1+ϵ22^M2 , (30)

with

 ^Mi=^Mi(Ψ)=−4π^ρ(Ψ)9λi . (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 Poisson-Laplace 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 second-order 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 non-monotonic dependence on concentration is revealed, with generally higher average eccentricities attained by low-concentration 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 mass-to-light 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

 ^Σ(^xP,^yP)=∫^zsp−^zsp^ρ(^rP)d^zP , (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 equally-spaced 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 semi-axes, 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 second-order 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 second-order 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, circular-averaged 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

 σ2(ψ)=25aγ(7/2,ψ)γ(5/2,ψ)=1a^σ2(ψ) , (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 one-to-one 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

 σ2P(^xP,^yP)=∫^zsp−^zspσ2(^rP)ρ(^rP)d^zPΣ(^xP,^yP)= 25a∫^zsp−^zspγ[7/2,ψ(^rP)]exp[ψ(^rP)]d^zP^Σ(^xP,^yP)=1a^σ2P(^xP,^yP) . (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 three-dimensional surface of these tidally truncated models. In § 2 of the present Paper, we have further emphasized the additional singularity that characterizes close-to-critical 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 (N-body, Fokker-Planck, Monte Carlo, gas), on the issue of the interplay between pressure anisotropy and tidal effects.

Since the first solutions of the Fokker-Planck equation by means of a Monte Carlo approach, as described by Hénon (1971) or Spitzer & Hart (1971), it has been shown that one-component isolated spherical clusters, starting from a variety of initial conditions (see Spitzer 1987; chap.4 and references therein), develop a core-halo 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 N-body simulations (e.g., see Giersz & Heggie 1997; Aarseth & Heggie 1998; Baumgardt & Makino 2003; Lee et al. 2006), anisotropic Fokker-Planck (e.g., see Takahashi et al. 1997), and Monte Carlo models (e.g., see Giersz 2001) of both single and multi-mass 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 time-dependent 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 Fokker-Planck 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 zeroth-order 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:

1. Two tidal regimes exist, namely of low and high-deformation, 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 self-consistent nonspherical initial conditions of numerical simulations aiming at studying in further details the evolution of configurations starting from sub-critical tidal equilibria.

2. 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” zero-velocity 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 low-concentration models. This statement agrees with a general trend noted by White & Shawl (1987) for the globular cluster system of our Galaxy.

3. 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.

4. 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.

5. 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 first-order models, reduce to simple functions of the coefficient . Average eccentricities have been calculated from the inertia tensor components, evaluated numerically.

6. 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.

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

8. 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 “extra-tidal” 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.

We wish to thank D. C. Heggie for several discussions and for his code, available within the software enviroment STARLAB by P. Hut, S. McMillan, J. Makino, and S. Portegies Zwart. We also thank D. Chernoff, L. Ciotti, and E. Vesperini for a number of useful comments and suggestions. Our code implementing the iteration method, described in Appendix A, makes use of the package S2kit 1.0 by P.J. Kostelec and D.N. Rockmore, and of FFTW 3.2.1 by M. Frigo and S. G. Johnson. The final part of this work was carried out at The Kavli Institute for Theoretical Physics, Santa Barbara, while we participated in the Program “Formation and Evolution of Globular Clusters”.

## 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 right-hand 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 second-order 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

 aΦ\tiny(ext)C(^r)=−94π^ρ(Ψ)∫Vd3^r′^ρ(^r′)|^r−^r′|=−9^ρ(Ψ)+∞∑l=012l+1+l∑m=−l^qlmYlm(θ,ϕ)^rl+1 , (B1)

where the multipole coefficients are defined as

 ^qlm=∫Vd3^r′Ylm(θ′,ϕ′)^r′l^ρ(^r′) , (B2)

with the integration to be performed in the volume with boundary surface defined by . This general expression can be compared with the second-order solution of the Laplace equation obtained in Paper I

 (aΦ\tiny(ext)C)(2)(^r) = [λ0+λ1ϵ+λ2ϵ22]1^r+[a20ϵ+b20ϵ22]Y20(θ,ϕ)^r3+[a22ϵ+b22ϵ22]Y22(θ,ϕ)^r3 (B3) +b40ϵ22Y40(θ,ϕ)^r5+b42ϵ22Y42(θ,ϕ)^r5+b44ϵ22Y44(θ,ϕ)^r5 ,

in order to determine the relation between the second-order multipole coefficients (i.e. calculated by means of the second-order expression for the density) and the matching constants that appear in Eq. (B3). Therefore, we find that

 λ0+λ1ϵ+λ2ϵ