Non-spherical quasi-relaxed stellar systems

# The construction of non-spherical models of quasi-relaxed stellar systems

G. Bertin and A.L. Varri Università degli Studi di Milano, Dipartimento di Fisica, via Celoria 16, I-20133 Milano, Italy; anna.varri@studenti.unimi.it
###### Abstract

Spherical models of collisionless but quasi-relaxed stellar systems have long been studied as a natural framework for the description of globular clusters. Here we consider the construction of self-consistent models under the same physical conditions, but including explicitly the ingredients that lead to departures from spherical symmetry. In particular, we focus on the effects of the tidal field associated with the hosting galaxy. We then take a stellar system on a circular orbit inside a galaxy represented as a “frozen” external field. The equilibrium distribution function is obtained from the one describing the spherical case by replacing the energy integral with the relevant Jacobi integral in the presence of the external tidal field. Then the construction of the model requires the investigation of a singular perturbation problem for an elliptic partial differential equation with a free boundary, for which we provide a method of solution to any desired order, with explicit solutions to two orders. We outline the relevant parameter space, thus opening the way to a systematic study of the properties of a two-parameter family of physically justified non-spherical models of quasi-relaxed stellar systems. The general method developed here can also be used to construct models for which the non-spherical shape is due to internal rotation. Eventually, the models will be a useful tool to investigate whether the shapes of globular clusters are primarily determined by internal rotation, by external tides, or by pressure anisotropy.

globular clusters: general — methods: analytical — stellar dynamics

## 1 Introduction

Large stellar systems can be studied as collisionless systems, by means of a one-star distribution function obeying the combined set of the collisionless Boltzmann equation and the Poisson equation, under the action of the self-consistent mean potential. For elliptical galaxies the relevant two-star relaxation times do actually exceed their age; an imprint of partial relaxation may be left at the time of their formation if we refer to a picture of formation via incomplete violent relaxation; Lynden-Bell 1967, van Albada 1982, but otherwise they should be thought of as truly collisionless systems, generally characterized by an anisotropic pressure tensor. In turn, for globular clusters the relevant relaxation times are typically shorter than their age, so that we may argue that for many of them the two-star relaxation processes have had enough time to act and to bring them close to a thermodynamically relaxed state, with their distribution function close to a Maxwellian. This line of arguments has led to the development of well-known dynamical models for globular clusters (King 1965; 1966).

King models are based on a quasi-Maxwellian isotropic distribution function in which a truncation prescription, continuous in phase space, is set heuristically to incorporate the presence of external tides; but otherwise, they are fully self-consistent (i.e., no external fields are actually considered) and perfectly spherical. Empirically, the simplification of spherical symmetry is encouraged by the fact that in general globular clusters have round appearance. Indeed, as a zero-th order description, these models have had remarkable success in applications to observed globular clusters (e.g., see Spitzer 1987; Djorgovski & Meylan 1994; McLaughlin & van der Marel 2005; and references therein). In recent years, great progress has been made in the acquisition of detailed quantitative information about the structure of these stellar systems, especially in relation to the measurement of the proper motions of thousands of individual stars (see van Leeuwen et al. 2000; McLaughlin et al. 2006), with the possibility of getting a direct 5-dimensional view of their phase space. Such progress calls for renewed efforts on the side of modeling. More general models would allow us to address the issue of the origin of the observed departures from spherical symmetry. In fact, it remains to be established which physical ingredient among rotation, pressure anisotropy, and tides is the primary cause of the flattening of globular clusters (e.g., see King 1961; Frenk & Fall 1982; Geyer et al. 1983; White & Shawl 1987; Davoust & Prugniel 1990; Han & Ryden 1994; Ryden 1996; Goodwin 1997; van den Bergh 2008; and references therein)

As in the case of the study of elliptical galaxies (e.g., see Bertin & Stiavelli 1993; and references therein), different approaches can be taken to the construction of models. Broadly speaking, two complementary paths can be followed. In the first, “descriptive” approach, under suitable geometrical (on the intrinsic shape) and dynamical (e.g., on the absence or presence of dark matter) hypotheses, the available data for an individual stellar system are imposed as constraints to derive the internal orbital structure (distribution function) most likely to correspond to the observations. This approach is often carried out in terms of codes that generalize a method introduced by Schwarzschild (1979); for an application to the globular cluster Cen, see van de Ven et al. (2006). In the second, “predictive” approach, one proposes a formation/evolution scenario in order to identify a physically justified distribution function for a wide class of objects, and then proceeds to investigate, by comparison with observations of several individual objects, whether the data support the general physical picture that has been proposed. Indeed, King models belong to this latter approach.

The purpose of this paper is to extend the description of quasi-relaxed stellar systems, so far basically limited to the spherical King models, to the non-spherical case. There are at least three different ways of extending spherical isotropic models of quasi-relaxed stellar systems (such as King models), by modifying the distribution function so as to include: (i) the explicit presence of a non-spherical tidal field; (ii) the presence of internal rotation; (iii) the presence of some pressure anisotropy. As noted earlier in this Introduction, these correspond to the physical ingredients that, separately, may be thought to be at the origin of the observed non-spherical shapes. We will thus focus on the construction of physically justified models, as an extension of the King models, in the presence of external tides and, briefly, on the extension of the models to the presence of rigid internal rotation. A first-order analysis of the triaxial tidal problem addressed in this paper was carried out by Heggie & Ramamani (1995) and the effect of a “frozen” tidal field on (initially) spherical King models was studied by Weinberg (1993) using N-body simulations. For the axisymmetric problem associated with the presence of internal rotation, a first-order analysis of a simply truncated Maxwellian distribution perturbed by a rigid rotation was given by Kormendy & Anand (1971); different models where differential rotation is considered were proposed by Prendergast & Tomer (1970) and by Wilson (1975), also in view of extensions to the presence of pressure anisotropy, which goes beyond the scope of this paper. Models that represent a direct generalization of the King family to the case with differential rotation have also been examined, with particular attention to their thermodynamic properties (Lagoute & Longaretti 1996; Longaretti & Lagoute 1996). In principle, the method of solution that we will present below can deal with the extension of other spherical isotropic models with finite size that are not of the King form e.g., see Woolley & Dickens 1962; Davoust 1977; see also the interesting suggestion by Madsen 1996.

The study of self-consistent collisionless equilibrium models has a long tradition not only in stellar dynamics, but also in plasma physics (e.g., see Harris 1962; Attico & Pegoraro 1999). We note that in both research areas a study in the presence of external fields, especially when the external field is bound to break the natural symmetry associated with the one-component problem, is only rarely considered.

The paper is organized as follows. Section 2 introduces the reference physical model, in which a globular cluster is imagined to move on a circular orbit inside a host galaxy treated as a frozen background field; the modified distribution function for such a cluster is then identified and the relevant parameter space defined. In Sect. 3 we set the mathematical problem associated with the construction of the related self-consistent models. For models generated by the spherical , Sect. 4 and Appendix A give the complete solution in terms of matched asymptotic expansions. Alternative methods of solutions are briefly discussed in Sect. 5. The concluding Section 6 gives a summary of the paper, with a short discussion of the results obtained. In Appendix B we show how the method developed in this paper can be applied to construct quasi-relaxed models flattened by rotation in the absence of external tides. In Appendix C we show how the method can be applied to other isotropic truncated models, different from King models.

Technically, the mathematical problem of a singular perturbation with a free boundary that is faced here is very similar to the problem noted in the theory of rotating stars, starting with Milne (see Tassoul 1978; Milne 1923; Chandrasekhar 1933; Krogdahl 1942; Chandrasekhar & Lebovitz 1962; Monaghan & Roxburgh 1965). The problem was initially dealt with inadequate tools; a satisfactory solution of the singular perturbation problem was obtained only later, by Smith (1975; 1976).

## 2 The physical model

### 2.1 The tidal potential

As a reference case, we consider an idealized model in which the center of mass of a globular cluster is imagined to move on a circular orbit of radius , characterized by orbital frequency , inside a host galaxy. For simplicity, we focus on the motion of the stars inside the globular cluster and model the galaxy, taken to have very large mass, by means of a frozen gravitational field (which we will call the galactic field, described by the potential ), with a given overall symmetry. This choice makes us ignore interesting effects that are generally present in the full interaction between a “satellite” and a galaxy; in a sense, we are taking a complementary view of an extremely complex dynamical situation, with respect to other investigations, such as those that lead to a discussion of the mechanism of dynamical friction in which the globular cluster or satellite is modeled as a rigid body and the stars of the galaxy are taken as the “live” component; see Bontekoe & van Albada 1987; Bertin et al. 2003; Arena & Bertin 2007; and references therein. Therefore, we will be initially following the picture of a restricted three-body problem, with one important difference, that the “secondary” is not treated as a point mass but as a “live” stellar system, described by the cluster mean-field potential . In this extremely simple orbital choice for the cluster center of mass, in the corotating frame the associated tidal field is time-independent and so we can proceed to the construction of a stationary dynamical model.

We consider the galactic potential to be spherically symmetric, i.e. , with , in terms of a standard set of Cartesian coordinates , so that . Let be the orbit plane of the center of mass of the cluster. We then introduce a local rotating frame of reference, so that the position vector is given by , with origin in the center of mass of the cluster and for which the x-axis points away from the center of the galactic field, the y-axis follows the direction of the cluster rotation in its orbit around the galaxy, and the z-axis is perpendicular to the orbit plane (according to the right-hand rule). In such rotating local frame, the relevant Lagrangian, describing the motion of a star belonging to the cluster, is cf. Chandrasekhar 1942, Eq. [5.510]:

 L=12{˙x2+˙y2+˙z2+Ω2[(R0+x)2+y2]+2Ω(R0+x)˙y−2Ω˙xy}−ΦG(R)−ΦC(x,y,z) , (1)

where and the terms responsible for centrifugal and Coriolis forces are explicitly displayed.

If we suppose that the size of the cluster is small compared to , we can adequately represent the galactic field by its linear approximation with respect to the local coordinates (the so-called “tidal approximation”). The corresponding equations of the motion for a single star in the rotating local frame are:

 ¨x−2Ω˙y−(4Ω2−κ2)x=−∂ΦC∂x , (2) ¨y+2Ω˙x=−∂ΦC∂y , (3) ¨z+Ω2z=−∂ΦC∂z , (4)

where is the epicyclic frequency at , given by . Note that the assumed symmetry for introduces a cancellation between the kinematic term and the gradient of the galactic potential and makes the vertical acceleration approximately equal to .

These equations admit an energy (isolating) integral of the motion, known as the Jacobi integral:

 H=12(˙x2+˙y2+˙z2)+ΦT+ΦC , (5)

where

 ΦT=12Ω2(z2−νx2) (6)

is the tidal potential. Here is a generally positive dimensionless coefficient.

Thus, at the level of single star orbits, we note that, in general, the tidal potential leads to a compression in the z-direction, a stretching in the x-direction, and leaves the y-direction untouched. The tidal potential is static, breaks the spherical symmetry, but is characterized by reflection symmetry with respect to the three natural coordinate planes; strictly speaking, the symmetry with respect to is applicable only in the limit of an infinitely massive host galaxy (see Spitzer 1987). In turn, we will see that the geometry of the tidal potential induces a non-spherical distortion of the cluster shape collectively, in particular an elongation along the x-axis and a compression along the z-axis. In practice, the numerical coefficient that defines quantitatively the induced distortion depends on the galactic potential. We recall that we have for a Keplerian potential, for a logarithmic potential, while for a Plummer model the dimensionless coefficient depends on the location of the circular orbit with respect to the model scale radius , with for a definition of the Plummer model see, e.g., Bertin 2000.

Different assumptions on the geometry of the galactic field can be treated with tools similar to those developed here, leading to a similar structure of the equations of the motion, with a slight modification of the tidal field. In particular, for an axisymmetric galactic field, the tidal potential differs from the one obtained here only by the z-term (Chandrasekhar 1942; Heggie & Hut 2003). This case is often considered, for example by referring to a globular cluster in circular orbit on the (axisymmetric) disk of our Galaxy (see Heggie & Ramamani 1995; Ernst et al. 2008), for which is then formulated in terms of the Oort constants.

In the physical model outlined in this Section, the typical dynamical time associated with the star orbits inside the cluster is much smaller than the (external) orbital time associated with . Therefore, in an asymptotic sense, the equilibrium configurations that we will construct in the rest of the paper can actually be generalized, with due qualifications, to more general orbits of the cluster inside a galaxy, provided we interpret the results that we are going to obtain as applicable only to a small piece of the cluster orbit.

### 2.2 The distribution function

As outlined in the Introduction, we wish to extend the description of quasi-relaxed stellar systems (so far basically limited to spherical models associated with distribution functions , dependent only on the single-star specific energy ) to the non spherical case, by including the presence of a non-spherical tidal field explicitly. Given the success of the spherical King models in the study of globular clusters, we will focus on the extension of models based on , which is defined as a “lowered” Maxwellian, continuous in phase space, with an energy cut-off which implies the existence of a boundary at the truncation radius .

Therefore, we will consider (partially) self-consistent models characterized by the distribution function:

 fK(H)=A[exp(−aH)−exp(−aH0)] (7)

if and otherwise, in terms of the Jacobi integral defined by Eq. (5). Here is the cut-off value for the Jacobi integral, while and are positive constants.

In velocity space, the inequality identifies a spherical region given by , where:

 ψ(r)=a{H0−[ΦC(r)+ΦT(x,z)]} (8)

is the dimensionless escape energy. Therefore, the boundary of the cluster is defined as the relevant zero-velocity surface by the condition and is given only implicitly by an equipotential (Hill) surface for the total potential ; in fact, its geometry depends on the properties of the tidal potential (of known characteristics; see Eq. [6]) and of the cluster potential (unknown a priori, to be determined as the solution of the associated Poisson equation).

The value of the cut-off potential should be chosen in such a way that the surface that defines the boundary is closed. The last (i.e., outermost) closed Hill surface is a critical surface, because it contains two saddle points that represent the Lagrangian points of the restricted three-body problem outlined in the previous subsection. From Eqs. (2)-(4), we see that such Lagrangian points are located symmetrically with respect to the origin of the local frame of reference and lie on the x-axis. Their distance from the origin is called the tidal radius, which we denote by , and can be determined from the condition:

 ∂ψ∂x(rT,0,0)=0 . (9)

If, as a zero-th order approximation, we use a simple Keplerian potential for the cluster potential , we recover the classical expression (e.g., see Spitzer 1987):

 r(0)T=(GMΩ2ν)1/3 , (10)

where is the total mass of the cluster.

As for the spherical King model, the density profile associated with the distribution function (7) is given by:

 ρ(ψ)=^Aeψγ(52,ψ)≡^A^ρ(ψ) , (11)

where . We recall that the incomplete gamma function has non-negative real value only in correspondence to a non-negative argument. In the following, we will denote the central density of the cluster by , where is the depth of the central potential well.

### 2.3 The parameter space

The models defined by are characterized by two physical scales (e.g., the two free constants and , or, correspondingly, the total mass and the central density of the cluster) and two dimensionless parameters. The first dimensionless parameter can be defined, as in the spherical King models, to measure the concentration of the cluster. We can thus consider the quantity , introduced at the end of the previous subsection, or we may refer to the commonly used concentration parameter:

 C=log(rtr/r0) , (12)

where is a scale length related to the size of the core and is the truncation radius of the spherical King model associated with the same value of the central potential well the relation between and is one-to-one; e.g., see Bertin 2000.

The second dimensionless parameter characterizes the strength of the (external) tidal field:

 ϵ≡Ω24πGρ0 . (13)

The definition arises naturally from the dimensionless formulation of the Poisson equation that describes the (partially) self-consistent problem (to be addressed in the next Section).

In principle, for a given choice of the dimensional scales ( and ) the truncation radius or the concentration parameter of a spherical King model can be set arbitrarily. In practice, the physical motivation of the models suggests that the truncation radius should be taken to be of the order of (and not exceed) the tidal radius , introduced in the previous subsection (see Eq. [9]). We may thus define an extension parameter, as the ratio between the truncation radius of the corresponding spherical model and the tidal radius :

 δ≡rtrrT . (14)

For a given value of the central potential well , there exists a (maximum) critical value for the tidal strength parameter, which we will denote by , corresponding to the maximum value for the extension parameter , which can be found by solving the system:

 ⎧⎪⎨⎪⎩∂ψ∂x(rT,0,0;ϵcr)=0ψ(rT,0,0;ϵcr)=0 . (15)

From this system, if we use the zero-th order Keplerian approximation for , we find that (see Spitzer 1987).

For our two-parameter family of models we thus expect two tidal regimes to exist. For models characterized by the pairs near the critical condition the tidal distortion should be maximal, while for models with pairs well below criticality only small departures from spherical symmetry should occur. A thorough exploration of the parameter space will be carried out in a separate paper (Varri & Bertin, in preparation). In closing, we note that the models proposed and studied by Heggie & Ramamani (1995) correspond to the pairs in parameter space that we have called critical.

## 3 The mathematical problem

The (partially) self-consistent models associated with the distribution function defined by Eq. (7) are constructed by solving the relevant Poisson equation. In terms of the dimensionless escape energy , given by Eq. (8), the Poisson equation (for ) can be written as:

 ∇2(ψ+aΦT)=−9r20ρρ0=−9r20^ρ(ψ)^ρ(Ψ) , (16)

where is the scale length introduced in Sect. 2.3 (see Eq. [12]). We then rescale the coordinates and introduce the dimensionless position vector , so that and , where we have made use of the tidal parameter introduced in Eq. (13). Therefore, the Poisson equation, for , can be written in dimensionless form as:

 ^∇2ψ=−9[^ρ(ψ)^ρ(Ψ)+ϵ(1−ν)] , (17)

while for negative values of we should refer to:

 ^∇2ψ=−9ϵ(1−ν) , (18)

i.e. the Laplace equation .

The mathematical problem is completed by specifying the appropriate boundary conditions. As for the spherical King models, we require regularity of the solution at the origin

 ψ(0)=Ψ , (19) ^∇ψ(0)=0 , (20)

 ψ+ϵT→aH0 , (21)

which corresponds to .

Poisson and Laplace domains are thus separated by the surface defined by which is unknown a priori; in other words, we have to solve an elliptic partial differential equation in a free boundary problem.

In the ordinary differential problem that characterizes the construction of spherical models with finite mass, the condition of vanishing cluster potential at large radii (together with the regularity conditions at the origin) overdetermines the problem, which can then be seen as an eigenvalue problem e.g., see Sect. 2.5 in Bertin & Stiavelli 1993. Indeed, for the King models the integration of the Poisson equation from the origin outwards, with “initial conditions” (19)-(20), sets the relation between the ratio and in order to meet the requirement (21), with thus playing the role of an “eigenvalue”.

In the more complex, three-dimensional situation that we are facing here, the existence of two different domains, internal (Poisson) and external (Laplace), suggests the use of the method of matched asymptotic expansions in order to obtain a uniform solution across the separation free surface. The solution in the internal and external domains are expressed as an asymptotic series with respect to the tidal parameter , which is assumed to be small (following the physical model described in the previous Section):

 ψ(int)(^r;ϵ)=∞∑k=01k!ψ(int)k(^r)ϵk , (22)
 ψ(ext)(^r;ϵ)=∞∑k=01k!ψ(ext)k(^r)ϵk , (23)

with spherical symmetry assumed for the zero-th order terms. The internal solution should obey the boundary conditions (19)-(20), while the external solution should satisfy Eq. (21). The two representations should be properly connected at the surface of the cluster.

On the other hand, for any small but finite value of the boundary, defined by , will be different from the unperturbed boundary, defined by , so that, for each of the two representations given above, there will be a small region in the vicinity of the surface of the cluster where the leading term is vanishingly small and actually smaller than the remaining terms of the formal series. Therefore, we expect the validity of the expansion to break down where the second term becomes comparable to the first, i.e where . This region can be considered as a boundary layer, which should be examined in “microscopic” detail by a suitable rescaling of the spatial coordinates and for which an adequate solution , expressed as a different asymptotic series, should be constructed. To obtain a uniformly valid solution over the entire space, an asymptotic matching is performed between the pairs (, ) and (, ), thus leading to a solution , obeying all the desired boundary conditions, in terms of three different, but matched, representations. This method of solution is basically the same method proposed by Smith (1975) for the analogous mathematical problem that arises in the determination of the structure of rigidly rotating fluid polytropes.

## 4 Solution in terms of matched asymptotic expansions

The complete solution to two significant orders in the tidal parameter is now presented. The formal solution to three orders is also displayed because of the requirements of the Van Dyke principle of asymptotic matching cf. Van Dyke 1975, Eq. [5.24] that we have adopted.

### 4.1 Internal region

If we insert the series (22) in the Poisson equation (17), under the conditions (19)-(20), we obtain an (infinite) set of Cauchy problems for . The problem for the zero-th order term (i.e., the unperturbed problem with ) is the one defining the construction of the spherical and fully self-consistent King models:

 ψ(int)0′′+2^rψ(int)0′=−9^ρ(ψ(int)0)^ρ(Ψ) , (24)

with and , where the symbol denotes derivative with respect to the argument . We recall that the truncation radius , which defines the boundary of the spherical models, is given implicitly by .

Let us introduce the quantities:

 Rj(^r;Ψ)≡9^ρ(Ψ)dj^ρdψj∣∣∣ψ(int)0 . (25)

These quantities depend on implicitly, through the function ; in turn, the dependence on is both explicit (through the term ) and implicit (because the function depends on the value of ). For convenience, we give the expression of the first terms of the sequence (cf. Eq. [11]): , , . Note that for , i.e. for , , , while for the quantities actually diverge. This is one more indication of the singular character of our perturbation analysis, which brings in some fractional power dependence on the perturbation parameter (see also expansion [46]).

Therefore, the equations governing the next two orders (for with ) can be written as:

 [^∇2+R1(^r;Ψ)]ψ(int)1=−9(1−ν) (26)
 [^∇2+R1(^r;Ψ)]ψ(int)2=−R2(^r;Ψ)(ψ(int)1)2 (27)

with and . The equation for is recorded in Appendix A.1, where we also describe the structure of the general equation for .

For any given order of the expansion, the operator acting on the function (see the left-hand side of Eqs. [26] and [27]) is the same, i.e. a Laplacian “shifted” by the function . If we thus expand every term in spherical harmonics111We use orthonormalized real spherical harmonics with Condon-Shortley phase; with respect to the toroidal angle , they are even for and odd otherwise.:

 ψ(int)k(^r)=∞∑l=0l∑m=−lψ(int)k,lm(^r)Ylm(θ,ϕ) , (28)

the three-dimensional differential problem is reduced to a one-dimensional (radial) problem, characterized by the following second order, linear ordinary differential operator:

 Dl=d2d^r2+2^rdd^r−l(l+1)^r2+R1(^r;Ψ) . (29)

In general, for a fixed value of , two independent solutions to the homogeneous problem are expected 222We note that , i.e. a numerical positive constant. Therefore, for the operator tends to the operator associated with the spherical Bessel functions of the first and of the second kind e.g., see Abramowitz & Stegun 1964, Eq. [10.1.1] for the equation and Eqs. [10.1.4] and [10.1.5] for the limiting values of the functions for small argument. to behave like and for . Because of the presence of , solutions to equations where appears have to be obtained numerically.

For (see Eq. [26]) we thus have to address the following problem. For , the relevant equation is:

 D0f00=−9(1−ν), (30)

where , with . Here we do not have to worry about including solutions to the associated homogeneous problem, because one of the two independent solutions would be singular at the origin and the other would be forced to vanish by the required condition at . For we have:

 Dlψ(int)1,lm=0 (31)

with . Both Eq. (31) and the associated boundary conditions are homogeneous. Therefore, the solution is undetermined by an -dependent multiplicative constant: , with for (the singular solution is excluded by the boundary conditions at the origin). Then the complete formal solution is:

 ψ(int)1(^r)=f00(^r)+∞∑l=1l∑m=−lAlmγl(^r)Ylm(θ,ϕ) , (32)

where the constants are ready to be determined by means of the asymptotic matching with at the boundary layer.

For (see Eq. [27]) the relevant equations are:

 Dlψ(int)2,lm=−R2(^r;Ψ)[ψ(int)12]lm , (33)

where on the right-hand side the function is that derived from the solution of the first order problem (which shows the progressive character of this method for the construction of solutions). In Appendix A.2 the equations for the six relevant harmonics are displayed explicitly. The boundary conditions to be imposed at the origin are again homogeneous: . For a fixed harmonic with , the general solution of Eq. (33) is the sum of a particular solution (which we will denote by ) and of a regular solution to the associated homogeneous problem given by Eq. (31) (which we will call , with the same functions introduced for the first order problem). Obviously, the particular solution exists only when Eq. (33) is non-homogeneous, i.e. only for those values of that correspond to a non-vanishing coefficient in the expansion of in spherical harmonics. As noted in the first order problem (), the associated homogeneous problem for has no non-trivial solution. Then we can express the complete solution as:

 ψ(int)2(^r)=g00(^r)+∞∑l=1l∑m=−l[glm(^r)+Blmγl(^r)]Ylm(θ,ϕ) , (34)

where are constants to be determined from the matching with the boundary layer.

Similarly, for the solution can be written as:

 ψ(int)3(^r)=h00(^r)+∞∑l=1l∑m=−l[hlm(^r)+Clmγl(^r)]Ylm(θ,ϕ) , (35)

where are particular solutions and are constants, again to be determined from the matching with the boundary layer.

Because the differential operator and the boundary conditions at the origin are the same for the reduced radial problem of every order, we have thus obtained the general structure of the solution for the internal region (see Appendix A.1).

### 4.2 External region

Here we first present the general solution and then proceed to set up the asymptotic series (23).

The solution to Eq. (18) describing the external region, i.e. in the Laplace domain, can be expressed as the sum of a particular solution () and of the solutions to the radial part of the Laplacian operator consistent with the boundary condition (21):

 ψ(ext)(^r)=α−λ^r−∞∑l=1l∑m=−lβlm^rl+1Ylm(θ,ϕ)−ϵT(^r) . (36)

Here we note that the tidal potential contributes only with spherical harmonics of order with even values of :

 T00(^r)=−3√π(ν−1)^r2 , (37) T20(^r)=3√π5(2+ν)^r2 , (38) T22(^r)=−3√3π5ν^r2 . (39)

At this point we can proceed to set up the asymptotic series, by expanding the constant coefficients , , and with respect to :

 α=aH0=α0+α1ϵ+12!α2ϵ2+.. , (40)
 λ=λ0+λ1ϵ+12!λ2ϵ2+.. , (41)
 βlm=almϵ+12!blmϵ2+13!clmϵ3+.. . (42)

The last expansion starts with a first order term because the density distribution of the unperturbed problem is spherically symmetric.

For convenience, we give the explicit expression of the external solution up to third order:

 ψ(ext)(^r)=α0−λ0^r+{α1−λ1^r−T00(^r)2√π−∞∑l=1l∑m=−l[alm^rl+1+Tlm(^r)]Ylm(θ,ϕ)}ϵ +12![α2−λ2^r−∞∑l=1l∑m=−lblm^rl+1Ylm(θ,ϕ)]ϵ2+13![α3−λ3^r−∞∑l=1l∑m=−lclm^rl+1Ylm(θ,ϕ)]ϵ3 . (43)

### 4.3 Boundary layer

The boundary layer is the region where the function becomes vanishingly small. Since the unperturbed gravitational field at the truncation radius is finite, , for any value of , based on a Taylor expansion of about we may argue that the region in which the series (22) breaks down can be defined by . In this boundary layer we thus introduce a suitable change of variables:

 η=^rtr−^rϵ , (44)

take the ordering , and thus rescale the solution by introducing the function . For positive values of the Poisson equation (17) thus becomes:

 ∂2τ∂η2−2ϵ^rtr−ϵη∂τ∂η+ϵ2(^rtr−ϵη)2Λ2τ=−9^ρ(Ψ)ϵ^ρ(ϵτ)−9ϵ2(1−ν) , (45)

where is the angular part of the Laplacian in spherical coordinates. For negative values of we can write a similar equation, corresponding to Eq. (18), which is obtained from Eq. (45) by dropping the term proportional to .

With the help of the asymptotic expansion for small argument of the incomplete gamma function e.g., see Bender& Orszag 1999, Eq. [6.2.5], we find:

 ^ρ(ϵτ)∼25τ5/2ϵ5/2+435τ7/2ϵ7/2+... , (46)

so that, within the boundary layer, the contribution of (which is the one that distinguishes the Poisson from the Laplace regime) becomes significant only beyond the tidal term, as a correction .

Therefore, up to we can write:

 τ=τ0+τ1ϵ+12!τ2ϵ2. (47)

To this order, which is required for a full solution up to of the global problem (see Eqs. [22] and [23]), by equating in Eq. (45) the first powers of separately, we obtain the relevant equations for the first three terms:

 ∂2τ0∂η2=0 , (48)
 ∂2τ1∂η2=2^rtr∂τ0∂η , (49)
 ∂2τ2∂η2=4^rtr[∂τ1∂η+η^rtr∂τ0∂η]−2^r2trΛ2τ0−18(1−ν) . (50)

The equations are easily integrated in the variable , to obtain the solutions:

 τ0=F0(θ,ϕ)η+G0(θ,ϕ) , (51)
 τ1=F0(θ,ϕ)^rtrη2+F1(θ,ϕ)η+G1(θ,ϕ) , (52)
 τ2=2F0(θ,ϕ)^r2trη3−13^r2trΛ2F0(θ,ϕ)η3+2F1(θ,ϕ)^rtrη2 −9(1−ν)η2−1^r2trΛ2G0(θ,ϕ)η2+F2(θ,ϕ)η+G2(θ,ϕ) . (53)

The six free angular functions that appear in the formal solutions will be determined by the matching procedure.

### 4.4 Asymptotic matching to two orders

In order to obtain the solution, we must perform separately the relevant matching for the pairs and . We follow the Van Dyke matching principle, which requires that we compare the second order expansion of the internal and external solutions with the third order expansion of the boundary layer solution. The full procedure is described in Appendix A.3.

To first order (i.e., up to in series [22] and [23]), from the matching of the pair we find the free angular functions of (51) and (52):

 F0(θ,ϕ)=−ψ(int)0′(^rtr) , (54)
 G0(θ,ϕ)=ψ(int)1(^rtr,θ,ϕ) , (55)
 F1(θ,ϕ)=−∂ψ(int)1∂^r(^rtr,θ,ϕ) , (56)
 G1(θ,ϕ)=12ψ(int)2(^rtr,θ,ϕ) . (57)

From the matching of the pair we connect to the same angular functions, thus proving that the matching to first order is equivalent to imposing continuity of the solution up to second order and of the first radial derivative up to first order. This allows us to determine the free constants that are present in the first two terms of (4.2) and in (32):

 α0=λ0^rtr , (58)
 λ0=^r2trψ(int)0′(^rtr) , (59)
 α1=f00(^rtr)+^rtrf′00(^rtr)+3T00(^rtr)2√π , (60)
 λ1=^r2trf′00(^rtr)+^rtrT00(^rtr)√π , (61)
 A2m=−5T2m(^rtr)^rtrγ′2(^rtr)+3γ2(^rtr) , (62)
 a2m=−^r3tr[A2mγ2(^rtr)+T2m(^rtr)] . (63)

Note that if , for every value of , and that the constants for are non-vanishing only for . The constants that identify the solution are thus expressed in terms of the values of the unperturbed field , of the “driving” tidal potential , and of the solutions and (see Eqs. [30] and [31]) taken at .

The boundary surface of the first order model is defined implicitly by , i.e. the spherical shape of the King model is modified by monopole and quadrupole contributions, which are even with respect to toroidal and poloidal angles and characterized by reflection symmetry with respect to the three natural coordinates planes. As might have been expected from the physical model, the spherical shape is thus modified only by spherical harmonics for which the tidal potential has non-vanishing coefficients. Mathematically, this is non-trivial, because the first order equation in the internal region Eq. (26) is non-homogeneous only for ; the quadrupole contribution to the internal solution is formally “hidden” by the use of the function (which includes the tidal potential) and is unveiled by the matching which demonstrates that with are non-vanishing.

The first order solution can be inserted into the right-hand side of Eq. (33) to generate non-homogeneous equations (and thus particular solutions) only for and corresponding positive and even values of (see Appendix A.2). We can thus proceed to construct the second order solution in the same way described above for the first order solution. From the matching of the pair we determine the missing angular functions:

 F2(θ,ϕ)=−∂ψ(int)2∂^r(^rtr,θ,ϕ) , (64)
 G2(θ,ϕ)=13ψ(int)3(^rtr,θ,ϕ) , (65)

which are then connected to the properties of by the matching of the pair . This is equivalent to imposing continuity of the solution up to third order and of the first radial derivative up to second order and leads to the determination of the free constants that appear in the third term of (4.2) and in (34):

 α2=g00(^rtr)+^rtrg′00(^rtr) , (66)
 λ2=^r2trg′00(^rtr) , (67)
 B2m=−^rtrg′2m(^rtr)+3g2m(^rtr)^rtrγ′2(^rtr)+3γ2(^rtr) , (68)
 b2m=−^r3tr[g2m(^rtr)+B2mγ2(^rtr)] , (69)
 B4m=−^rtrg′4m(^rtr)+5g4m(^rtr)^rtrγ′4(^rtr)+5γ4(^rtr) , (70)
 b4m=−^r5tr[g4m(^rtr)+B4mγ4(^rtr)] . (71)

Here if for every value of ; the only non-vanishing constants with are those with even .

Therefore, the second order solution has non-vanishing contributions only for , i.e. for those harmonics for which the particular solution to Eq. (27) is non-trivial. By induction, it can be proved (see Appendix A.4) that the -th order solution is characterized by harmonics with corrisponding positive and even values of . In reality, the discussion of the matching to higher orders () would require a re-definition of the boundary layer, because the density contribution on the right-hand side of Eq. (45) (for positive values of ) comes into play. The asymptotic matching procedure carries through also in this more complex case but, for simplicity, is omitted here. We should also keep in mind that in an asymptotic analysis the inclusion of higher order terms does not necessarily lead to better accuracy in the solution; the optimal truncation in the asymptotic series depends on the value of the expansion parameter (in this case, on the value of ) and has to be judged empirically.

In conclusion, starting from a given value of the King concentration parameter and from a given strength of the tidal field , the uniform triaxial solution is constructed by numerically integrating Eqs. (24), (30), (31), and (33) and by applying the constants derived in this subsection to the asymptotic series expansion (22)-(23). The numerical integrations can be performed efficiently by means of standard Runge-Kutta routines. The boundary surface of the model is thus defined by , while the internal density distribution is given by , with the function defined by Eq. (11). Any other “observable” quantity can be reconstructed by suitable integration in phase space of the distribution function defined by Eq. (7), with defined by Eq. (5), and .

In Fig. 1 and Fig. 2 we illustrate the main characteristics of one triaxial model constructed with the method described in this Section.

## 5 Alternative methods of solution

### 5.1 The method of strained coordinates

The mathematical problem described in Sect. 3 can also be solved by the method of strained coordinates, an alternative method usually applied to non-linear hyperbolic differential equations e.g., see Van Dyke 1975, Chapter 6 and considered by Smith (1976) in the solution of the singular free-boundary perturbation problem that arises in the study of rotating polytropes.

Starting from a series representation of the form (22) and (23) for the solution defined in the Poisson and Laplace domains, respectively, a transformation is considered from spherical coordinates to “strained coordinates” :

 ^r=s+ϵ^r1(s,p,q)+12ϵ2^r2(s,p,q)+... θ=p ϕ=q , (72)

where are initially unspecified straining functions. We note that the zero-th order problem is defined by the same Eq. (24) with the same boundary conditions but with the variable replaced by . The unperturbed spherical boundary in the strained space is defined by , where . To each order, the effective boundary of the perturbed configuration remains described by the surface , while in physical coordinates the truncation radius actually changes as a result of the straining functions that are determined progressively.

The Laplacian expressed in the new coordinates, , can be written as an asymptotic series: