text:
# General solution of the Jeans equations for triaxial galaxies with separable potentials

## 1 Introduction

Much has been learned about the mass distribution and internal dynamics of galaxies by modeling their observed kinematics with solutions of the Jeans equations (e.g., Binney & Tremaine 1987). These are obtained by taking velocity moments of the collisionless Boltzmann equation for the phase-space distribution function $`f`$, and connect the second moments (or the velocity dispersions, if the mean streaming motion is known) directly to the density and the gravitational potential of the galaxy, without the need to know $`f`$. In nearly all cases there are fewer Jeans equations than velocity moments, so that additional assumptions have to be made about the degree of anisotropy. Furthermore, the resulting second moments may not correspond to a physical distribution function $`f0`$. These significant drawbacks have not prevented wide application of the Jeans approach to the kinematics of galaxies, even though the results need to be interpreted with care. Fortunately, efficient analytic and numerical methods have been developed in the past decade to calculate the full range of distribution functions $`f`$ that correspond to spherical or axisymmetric galaxies, and to fit them directly to kinematic measurements (e.g., Gerhard 1993; Qian et al. 1995; Rix et al. 1997; van der Marel et al. 1998). This has provided, for example, accurate intrinsic shapes, mass-to-light ratios, and central black hole masses (e.g., Verolme et al. 2002; Gebhardt et al. 2003).

Many galaxy components are not spherical or axisymmetric, but have triaxial shapes (Binney 1976, 1978). These include early-type bulges, bars, and giant elliptical galaxies. In this geometry, there are three Jeans equations, but little use has been made of them, as they contain six independent second moments, three of which have to be chosen ad-hoc (see, e.g., Evans, Carollo & de Zeeuw 2000). At the same time, not much is known about the range of physical solutions, as very few distribution functions have been computed, and even fewer have been compared with kinematic data (but see Zhao 1996).

An exception is provided by the special set of triaxial mass models that have a gravitational potential of Stäckel form. In these systems, the Hamilton–Jacobi equation separates in orthogonal curvilinear coordinates (Stäckel 1891), so that all orbits have three exact integrals of motion, which are quadratic in the velocities. The associated mass distributions can have arbitrary central axis ratios and a large range of density profiles, but they all have cores rather than central density cusps, which implies that they do not provide perfect fits to galaxies (de Zeeuw, Peletier & Franx 1986). Even so, they capture much of the rich internal dynamics of large elliptical galaxies (de Zeeuw 1985a, hereafter Z85; Statler 1987, 1991; Arnold, de Zeeuw & Hunter 1994). Numerical and analytic distribution functions have been constructed for these models (e.g., Bishop 1986; Statler 1987; Dejonghe & de Zeeuw 1988; Hunter & de Zeeuw 1992, hereafter HZ92; Mathieu & Dejonghe 1999), and their projected properties have been used to provide constraints on the intrinsic shapes of individual galaxies (e.g., Statler 1994a, b; Statler & Fry 1994; Statler, DeJonghe & Smecker-Hane 1999; Bak & Statler 2000; Statler 2001).

The Jeans equations for triaxial Stäckel systems have received little attention. This is remarkable, as Eddington (1915) already knew that the velocity ellipsoid in these models is everywhere aligned with the confocal ellipsoidal coordinate system in which the motion separates. This means that there are only three, and not six, non-vanishing second-order velocity moments in these coordinates, so that the Jeans equations form a closed system. However, Eddington, and later Chandrasekhar (1939, 1940), did not study the velocity moments, but instead assumed a form for the distribution function, and then determined which potentials are consistent with it. Lynden–Bell (1960) was the first to derive the explicit form of the Jeans equations for the triaxial Stäckel models. He showed that they constitute a highly symmetric set of three first-order partial differential equations (PDEs) for three unknowns, each of which is a function of the three confocal ellipsoidal coordinates, but he did not derive solutions. When it was realized that the orbital structure in the triaxial Stäckel models is very similar to that in the early numerical models for triaxial galaxies with cores (Schwarzschild 1979; Z85), interest in the second moments increased, and the Jeans equations were solved for a number of special cases. These include the axisymmetric limits and elliptic discs (Dejonghe & de Zeeuw 1988; Evans & Lynden–Bell 1989, hereafter EL89), triaxial galaxies with only thin tube orbits (HZ92), and, most recently, the scale-free limit (Evans et al. 2000). In all these cases the equations simplify to a two-dimensional problem, which can be solved with standard techniques after recasting two first-order equations into a single second-order equation in one dependent variable. However, these techniques do not carry over to a single third-order equation in one dependent variable, which is the best that one could expect to have in the general case. As a result, the general case has remained unsolved.

In this paper, we first present an alternative solution method for the two-dimensional limiting cases which does not use the standard approach, but instead uses superpositions of singular solutions. We show that this approach can be extended to three dimensions, and provides the general solution for the triaxial case in closed form, which we give explicitly. We will apply our solutions in a follow-up paper, and will use them together with the mean streaming motions (Statler 1994a) to study the properties of the observed velocity and dispersion fields of triaxial galaxies.

In §2, we define our notation and derive the Jeans equations for the triaxial Stäckel models in confocal ellipsoidal coordinates, together with the continuity conditions. We summarise the limiting cases, and show that the Jeans equations for all the cases with two degrees of freedom correspond to the same two-dimensional problem. We solve this problem in §3, first by employing a standard approach with a Riemann–Green function, and then via the singular solution superposition method. We also discuss the choice of boundary conditions in detail. We relate our solution to that derived by EL89 in Appendix A, and explain why it is different. In §4, we extend the singular solution approach to the three-dimensional problem, and derive the general solution of the Jeans equations for the triaxial case. It contains complete (hyper)elliptic integrals, which we express as single quadratures that can be numerically evaluated in a straightforward way. We summarise our conclusions in §5.

## 2 The Jeans equations for separable models

We first summarise the essential properties of the triaxial Stäckel models in confocal ellipsoidal coordinates. Further details can be found in Z85. We show that for these models the mixed second-order velocity moments vanish, so that the Jeans equations form a closed system. We derive the Jeans equations and find the corresponding continuity conditions for the general case of a triaxial galaxy. We then give an overview of the limiting cases and show that solving the Jeans equations for the various cases with two degrees of freedom reduces to an equivalent two-dimensional problem.

### 2.1 Triaxial Stäckel models

We define confocal ellipsoidal coordinates ($`\lambda ,\mu ,\nu `$) as the three roots for $`\tau `$ of

$$\frac{x^2}{\tau +\alpha }+\frac{y^2}{\tau +\beta }+\frac{z^2}{\tau +\gamma }=1,$$
(2.1)
with ($`x,y,z`$) the usual Cartesian coordinates, and with constants $`\alpha ,\beta `$ and $`\gamma `$ such that $`\gamma \nu \beta \mu \alpha \lambda `$. For each point ($`x,y,z`$), there is a unique set ($`\lambda ,\mu ,\nu `$), but a given combination ($`\lambda ,\mu ,\nu `$) generally corresponds to eight different points ($`\pm x,\pm y,\pm z`$). We assume all three-dimensional Stäckel models in this paper to be likewise eightfold symmetric.

Surfaces of constant $`\lambda `$ are ellipsoids, and surfaces of constant $`\mu `$ and $`\nu `$ are hyperboloids of one and two sheets, respectively (Fig. 1). The confocal ellipsoidal coordinates are approximately Cartesian near the origin and become a conical coordinate system at large radii with a system of spheres together with elliptic and hyperbolic cones (Fig. 3). At each point, the three coordinate surfaces are perpendicular to each other. Therefore, the line element is of the form $`ds^2=P^2d\lambda ^2+Q^2d\mu ^2+R^2d\nu ^2`$, with the metric coefficients

$`P^2`$ $`=`$ $`{\displaystyle \frac{(\lambda \mu )(\lambda \nu )}{4(\lambda +\alpha )(\lambda +\beta )(\lambda +\gamma )}},`$
$`Q^2`$ $`=`$ $`{\displaystyle \frac{(\mu \nu )(\mu \lambda )}{4(\mu +\alpha )(\mu +\beta )(\mu +\gamma )}},`$ (2.2)
$`R^2`$ $`=`$ $`{\displaystyle \frac{(\nu \lambda )(\nu \mu )}{4(\nu +\alpha )(\nu +\beta )(\nu +\gamma )}}.`$
We restrict attention to models with a gravitational potential $`V_S(\lambda ,\mu ,\nu )`$ of Stäckel form (Weinacht 1924)

$$V_S=\frac{F(\lambda )}{(\lambda \mu )(\lambda \nu )}\frac{F(\mu )}{(\mu \nu )(\mu \lambda )}\frac{F(\nu )}{(\nu \lambda )(\nu \mu )},$$
(2.3)
where $`F(\tau )`$ is an arbitrary smooth function.

Adding any linear function of $`\tau `$ to $`F(\tau )`$ changes $`V_S`$ by at most a constant, and hence has no effect on the dynamics. Following Z85, we use this freedom to write

$$F(\tau )=(\tau +\alpha )(\tau +\gamma )G(\tau ),$$
(2.4)
where $`G(\tau )`$ is smooth. It equals the potential along the intermediate axis. This choice will simplify the analysis of the large radii behaviour of the various limiting cases.111Other, equivalent, choices include $`F(\tau )=(\tau +\alpha )(\tau +\gamma )G(\tau )`$ by HZ92, and $`F(\tau )=(\tau +\alpha )(\tau +\beta )U(\tau )`$ by de Zeeuw et al. (1986), with $`U(\tau )`$ the potential along the short axis.

The density $`\rho _S`$ that corresponds to $`V_S`$ can be found from Poisson’s equation or by application of Kuzmin’s (1973) formula (see de Zeeuw 1985b). This formula shows that, once we have chosen the central axis ratios and the density along the short axis, the mass model is fixed everywhere by the requirement of separability. For centrally concentrated mass models, $`V_S`$ has the $`x`$-axis as long axis and the $`z`$-axis as short axis. In most cases this is also true for the associated density (de Zeeuw et al. 1986).

### 2.2 Velocity moments

A stellar system is completely described by its distribution function (DF), which in general is a time-dependent function $`f`$ of the six phase-space coordinates ($`𝐱,𝐯`$). Assuming the system to be in equilibrium ($`df/dt=0`$) and in steady-state ($`f/t=0`$), the DF is independent of time $`t`$ and satisfies the (stationary) collisionless Boltzmann equation (CBE). Integration of the DF over all velocities yields the zeroth-order velocity moment, which is the density $`\rho `$ of the stellar system. The first- and second-order velocity moments are defined as

$`v_i(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle v_if(𝐱,𝐯)\mathrm{d}^3v},`$
$`v_iv_j(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle v_iv_jf(𝐱,𝐯)\mathrm{d}^3v},`$
where $`i,j=1,2,3`$. The streaming motions $`v_i`$ together with the symmetric second-order velocity moments $`v_iv_j`$ provide the velocity dispersions $`\sigma _{ij}^2=v_iv_jv_iv_j`$.

The continuity equation that results from integrating the CBE over all velocities, relates the streaming motion to the density $`\rho `$ of the system. Integrating the CBE over all velocities after multiplication by each of the three velocity components, provides the Jeans equations, which relate the second-order velocity moments to $`\rho `$ and $`V`$, the potential of the system. Therefore, if the density and potential are known, we in general have one continuity equation with three unknown first-order velocity moments and three Jeans equations with six unknown second-order velocity moments.

The potential (2.3) is the most general form for which the Hamilton–Jacobi equation separates (Stäckel 1890; Lynden–Bell 1962b; Goldstein 1980). All orbits have three exact isolating integrals of motion, which are quadratic in the velocities (e.g., Z85). It follows that there are no irregular orbits, so that Jeans’ (1915) theorem is strictly valid (Lynden–Bell 1962a; Binney 1982) and the DF is a function of the three integrals. The orbital motion is a combination of three independent one-dimensional motions — either an oscillation or a rotation — in each of the three ellipsoidal coordinates. Different combinations of rotations and oscillations result in four families of orbits in triaxial Stäckel models (Kuzmin 1973; Z85): inner (I) and outer (O) long-axis tubes, short (S) axis tubes and box orbits. Stars on box orbits carry out an oscillation in all three coordinates, so that they provide no net contribution to the mean streaming. Stars on I- and O-tubes carry out a rotation in $`\nu `$ and those on S-tubes a rotation in $`\mu `$, and oscillations in the other two coordinates. The fractions of clockwise and counterclockwise stars on these orbits may be unequal. This means that each of the tube families can have at most one nonzero first-order velocity moment, related to $`\rho `$ by the continuity equation. Statler (1994a) used this property to construct velocity fields for triaxial Stäckel models. It is not difficult to show by similar arguments (e.g., HZ92) that all mixed second-order velocity moments also vanish

$$v_\lambda v_\mu =v_\mu v_\nu =v_\nu v_\lambda =0.$$
(2.6)
Eddington (1915) already knew that in a potential of the form (2.3), the axes of the velocity ellipsoid at any given point are perpendicular to the coordinate surfaces, so that the mixed second-order velocity moments are zero. We are left with three second-order velocity moments, $`v_\lambda ^2`$, $`v_\mu ^2`$ and $`v_\nu ^2`$, related by three Jeans equations.

### 2.3 The Jeans equations

The Jeans equations for triaxial Stäckel models in confocal ellipsoidal coordinates were first derived by Lynden–Bell (1960). We give an alternative derivation here, using the Hamilton equations.

We first write the DF as a function of ($`\lambda ,\mu ,\nu `$) and the conjugate momenta

$$p_\lambda =P^2\frac{d\lambda }{dt},p_\mu =Q^2\frac{d\mu }{dt},p_\nu =R^2\frac{d\nu }{dt},$$
(2.7)
with the metric coefficients $`P`$, $`Q`$ and $`R`$ given in (2.1). In these phase-space coordinates the steady-state CBE reads

$$\frac{d\tau }{dt}\frac{f}{\tau }+\frac{dp_\tau }{dt}\frac{f}{p_\tau }=0,$$
(2.8)
where we have used the summation convention with respect to $`\tau =\lambda ,\mu ,\nu `$. The Hamilton equations are

$$\frac{d\tau }{dt}=\frac{H}{p_\tau },\frac{dp_\tau }{dt}=\frac{H}{\tau },$$
(2.9)
with the Hamiltonian defined as

$$H=\frac{p_\lambda ^2}{2P^2}+\frac{p_\mu ^2}{2Q^2}+\frac{p_\nu ^2}{2R^2}+V(\lambda ,\mu ,\nu ).$$
(2.10)
The first Hamilton equation in (2.9) defines the momenta (2.7) and gives no new information. The second gives

$$\frac{dp_\lambda }{dt}=\frac{p_\lambda ^2}{P^3}\frac{P}{\lambda }+\frac{p_\mu ^2}{Q^3}\frac{Q}{\lambda }+\frac{p_\nu ^2}{R^3}\frac{R}{\lambda }\frac{V}{\lambda },$$
(2.11)
and similar for $`p_\mu `$ and $`p_\nu `$ by replacing the derivatives with respect to $`\lambda `$ by derivatives to $`\mu `$ and $`\nu `$, respectively.

We assume the potential to be of the form $`V_S`$ defined in (2.3), and we substitute (2.7) and (2.11) in the CBE (2.8). We multiply this equation by $`p_\lambda `$ and integrate over all momenta. The mixed second moments vanish (2.6), so that we are left with

$$\begin{array}{c}\frac{3fp_\lambda ^2}{P^3}\frac{P}{\lambda }+\frac{fp_\mu ^2}{Q^3}\frac{Q}{\lambda }+\frac{fp_\nu ^2}{R^3}\frac{R}{\lambda }\hfill \\ \hfill \frac{1}{P^2}\frac{}{\lambda }fp_\lambda ^2f\frac{V_S}{\lambda }=0,\end{array}$$
(2.12)
where we have defined the moments

$`f`$ $``$ $`{\displaystyle f\mathrm{d}^3p}=PQR\rho ,`$
$`fp_\lambda ^2`$ $``$ $`{\displaystyle p_\lambda ^2f\mathrm{d}^3p}=P^3QRT_{\lambda \lambda },`$
with the diagonal components of the stress tensor

$$T_{\tau \tau }(\lambda ,\mu ,\nu )\rho v_\tau ^2,\tau =\lambda ,\mu ,\nu .$$
(2.14)
The moments $`fp_\mu ^2`$ and $`fp_\nu ^2`$ follow from $`fp_\lambda ^2`$ by cyclic permutation $`\lambda \mu \nu \lambda `$, for which $`PQRP`$. We substitute the definitions (2.3) in eq. (2.12) and carry out the partial differentiation in the fourth term. The first term in (2.12) then cancels, and, after rearranging the remaining terms and dividing by $`PQR`$, we obtain

$$\frac{T_{\lambda \lambda }}{\lambda }+\frac{T_{\lambda \lambda }T_{\mu \mu }}{Q}\frac{Q}{\lambda }+\frac{T_{\lambda \lambda }T_{\nu \nu }}{R}\frac{R}{\lambda }=\rho \frac{V_S}{\lambda }.$$
(2.15)
Substituting the metric coefficients (2.1) and carrying out the partial differentiations results in the Jeans equations

$$\frac{T_{\lambda \lambda }}{\lambda }+\frac{T_{\lambda \lambda }T_{\mu \mu }}{2(\lambda \mu )}+\frac{T_{\lambda \lambda }T_{\nu \nu }}{2(\lambda \nu )}=\rho \frac{V_S}{\lambda },$$
(2.16a)
$$\frac{T_{\mu \mu }}{\mu }+\frac{T_{\mu \mu }T_{\nu \nu }}{2(\mu \nu )}+\frac{T_{\mu \mu }T_{\lambda \lambda }}{2(\mu \lambda )}=\rho \frac{V_S}{\mu },$$
(2.16b)
$$\frac{T_{\nu \nu }}{\nu }+\frac{T_{\nu \nu }T_{\lambda \lambda }}{2(\nu \lambda )}+\frac{T_{\nu \nu }T_{\mu \mu }}{2(\nu \mu )}=\rho \frac{V_S}{\nu },$$
(2.16c)
where the equations for $`\mu `$ and $`\nu `$ follow from the one for $`\lambda `$ by cyclic permutation. These equations are identical to those derived by Lynden–Bell (1960).

In self-consistent models, the density $`\rho `$ must equal $`\rho _S`$, with $`\rho _S`$ related to the potential $`V_S`$ (2.3) by Poisson’s equation. The Jeans equations, however, do not require self-consistency. Hence, we make no assumptions on the form of the density other than that it is triaxial, i.e., a function of $`(\lambda ,\mu ,\nu )`$, and that it tends to zero at infinity. The resulting solutions for the stresses $`T_{\tau \tau }`$ do not all correspond to physical distribution functions $`f0`$. The requirement that the $`T_{\tau \tau }`$ are non-negative removes many (but not all) of the unphysical solutions.

### 2.4 Continuity conditions

We saw in §2.2 that the velocity ellipsoid is everywhere aligned with the confocal ellipsoidal coordinates. When $`\lambda \alpha `$, the ellipsoidal coordinate surface degenerates into the area inside the focal ellipse (Fig. 2). The area outside the focal ellipse is labeled by $`\mu =\alpha `$. Hence, $`T_{\lambda \lambda }`$ is perpendicular to the surface inside and $`T_{\mu \mu }`$ is perpendicular to the surface outside the focal ellipse. On the focal ellipse, i.e. when $`\lambda =\mu =\alpha `$, both stress components therefore have to be equal. Similarly, $`T_{\mu \mu }`$ and $`T_{\nu \nu }`$ are perpendicular to the area inside ($`\mu =\beta `$) and outside ($`\nu =\beta `$) the two branches of the focal hyperbola, respectively, and have to be equal on the focal hyperbola itself ($`\mu =\nu =\beta `$). This results in the following two continuity conditions

$$T_{\lambda \lambda }(\alpha ,\alpha ,\nu )=T_{\mu \mu }(\alpha ,\alpha ,\nu ),$$
(2.17a)
$$T_{\mu \mu }(\lambda ,\beta ,\beta )=T_{\nu \nu }(\lambda ,\beta ,\beta ).$$
(2.17b)
These conditions not only follow from geometrical arguments, but are also precisely the conditions necessary to avoid singularities in the Jeans equations (2.16) when $`\lambda =\mu =\alpha `$ and $`\mu =\nu =\beta `$. For the sake of physical understanding, we will also obtain the corresponding continuity conditions by geometrical arguments for the limiting cases that follow.

### 2.5 Limiting cases

When two or all three of the constants $`\alpha `$, $`\beta `$ or $`\gamma `$ are equal, the triaxial Stäckel models reduce to limiting cases with more symmetry and thus with fewer degrees of freedom. We show in §2.6 that solving the Jeans equations for all the models with two degrees of freedom reduces to the same two-dimensional problem. EL89 first solved this generalised problem and applied it to the disc, oblate and prolate case. Evans et al. (2000) showed that the large radii case with scale-free DF reduces to the problem solved by EL89. We solve the same problem in a different way in §3, and obtain a simpler expression than EL89. In order to make application of the resulting solution straightforward, and to define a unified notation, we first give an overview of the limiting cases.

#### 2.5.1 Oblate spheroidal coordinates: prolate potentials

When $`\gamma =\beta `$, the coordinate surfaces for constant $`\lambda `$ and $`\mu `$ reduce to oblate spheroids and hyperboloids of revolution around the $`x`$-axis. Since the range of $`\nu `$ is zero, it cannot be used as a coordinate. The hyperboloids of two sheets are now planes containing the $`x`$-axis. We label these planes by an azimuthal angle $`\chi `$, defined as $`\mathrm{tan}\chi =z/y`$. In these oblate spheroidal coordinates ($`\lambda ,\mu ,\chi `$) the potential $`V_S`$ has the form (cf. Lynden–Bell 1962b)

$$V_S=\frac{f(\lambda )f(\mu )}{\lambda \mu }\frac{g(\chi )}{(\lambda +\beta )(\mu +\beta )},$$
(2.18)
where the function $`g(\chi )`$ is arbitrary, and $`f(\tau )=(\tau +\alpha )G(\tau )`$, with $`G(\tau )`$ as in eq. (2.4). The denominator of the second term is proportional to $`y^2+z^2`$, so that these potentials are singular along the entire $`x`$-axis unless $`g(\chi )0`$. In this case, the potential is prolate axisymmetric, and the associated density $`\rho _S`$ is generally prolate as well (de Zeeuw et al. 1986).

The Jeans equations (2.16) reduce to

$`{\displaystyle \frac{T_{\lambda \lambda }}{\lambda }}+{\displaystyle \frac{T_{\lambda \lambda }T_{\mu \mu }}{2(\lambda \mu )}}+{\displaystyle \frac{T_{\lambda \lambda }T_{\chi \chi }}{2(\lambda +\beta )}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\lambda }},`$
$`{\displaystyle \frac{T_{\mu \mu }}{\mu }}+{\displaystyle \frac{T_{\mu \mu }T_{\lambda \lambda }}{2(\mu \lambda )}}+{\displaystyle \frac{T_{\mu \mu }T_{\chi \chi }}{2(\mu +\beta )}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\mu }},`$ (2.19)
$`{\displaystyle \frac{T_{\chi \chi }}{\chi }}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\chi }}.`$
The continuity condition (2.17a) still holds, except that the focal ellipse has become a focal circle. For $`\mu =\beta `$, the one-sheeted hyperboloid degenerates into the $`x`$-axis, so that $`T_{\mu \mu }`$ is perpendicular to the $`x`$-axis and coincides with $`T_{\chi \chi }`$. This gives the following two continuity conditions

$`T_{\lambda \lambda }(\alpha ,\alpha ,\chi )`$ $`=`$ $`T_{\mu \mu }(\alpha ,\alpha ,\chi ),`$
$`T_{\mu \mu }(\lambda ,\beta ,\chi )`$ $`=`$ $`T_{\chi \chi }(\lambda ,\beta ,\chi ).`$
By integrating along characteristics, Hunter et al. (1990) obtained the solution of (2.5.1) for the special prolate models in which only the thin I- and O-tube orbits are populated, so that $`T_{\mu \mu }0`$ and $`T_{\lambda \lambda }0`$, respectively (cf. §2.5.6).

#### 2.5.2 Prolate spheroidal coordinates: oblate potentials

When $`\beta =\alpha `$, we cannot use $`\mu `$ as a coordinate and replace it by the azimuthal angle $`\varphi `$, defined as $`\mathrm{tan}\varphi =y/x`$. Surfaces of constant $`\lambda `$ and $`\nu `$ are confocal prolate spheroids and two-sheeted hyperboloids of revolution around the $`z`$-axis. The prolate spheroidal coordinates ($`\lambda ,\varphi ,\nu `$) follow from the oblate spheroidal coordinates ($`\lambda ,\mu ,\chi `$) by taking $`\mu \nu `$, $`\chi \varphi `$ and $`\beta \alpha \gamma `$. The potential $`V_S(\lambda ,\varphi ,\nu )`$ is (cf. Lynden–Bell 1962b)

$$V_S=\frac{f(\lambda )f(\nu )}{\lambda \nu }\frac{g(\varphi )}{(\lambda +\alpha )(\nu +\alpha )}.$$
(2.21)
In this case, the denominator of the second term is proportional to $`R^2=x^2+y^2`$, so that the potential is singular along the entire $`z`$-axis, unless $`g(\varphi )`$ vanishes. When $`g(\varphi )0`$, the potential is oblate, and the same is generally true for the associated density $`\rho _S`$.

The Jeans equations (2.16) reduce to

$`{\displaystyle \frac{T_{\lambda \lambda }}{\lambda }}+{\displaystyle \frac{T_{\lambda \lambda }T_{\varphi \varphi }}{2(\lambda +\alpha )}}+{\displaystyle \frac{T_{\lambda \lambda }T_{\nu \nu }}{2(\lambda \nu )}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\lambda }},`$
$`{\displaystyle \frac{T_{\varphi \varphi }}{\varphi }}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\varphi }}.`$ (2.22)
$`{\displaystyle \frac{T_{\nu \nu }}{\nu }}+{\displaystyle \frac{T_{\nu \nu }T_{\lambda \lambda }}{2(\nu \lambda )}}+{\displaystyle \frac{T_{\nu \nu }T_{\varphi \varphi }}{2(\nu +\alpha )}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\nu }}.`$
For $`\lambda =\alpha `$, the prolate spheroidal coordinate surfaces reduce to the part of the $`z`$-axis between the foci. The part beyond the foci is reached if $`\nu =\alpha `$. Hence, in this case, $`T_{\lambda \lambda }`$ is perpendicular to part of the $`z`$-axis between, and $`T_{\nu \nu }`$ is perpendicular to the part of the $`z`$-axis beyond the foci. They coincide at the foci ($`\lambda =\nu =\alpha `$), resulting in one continuity condition. Two more follow from the fact that $`T_{\varphi \varphi }`$ is perpendicular to the (complete) $`z`$-axis, and thus coincides with $`T_{\lambda \lambda }`$ and $`T_{\nu \nu }`$ on the part between and beyond the foci, respectively:

$`T_{\lambda \lambda }(\alpha ,\varphi ,\alpha )`$ $`=`$ $`T_{\nu \nu }(\alpha ,\varphi ,\alpha ),`$
$`T_{\lambda \lambda }(\alpha ,\varphi ,\nu )`$ $`=`$ $`T_{\varphi \varphi }(\alpha ,\varphi ,\nu ),`$ (2.23)
$`T_{\nu \nu }(\lambda ,\varphi ,\alpha )`$ $`=`$ $`T_{\varphi \varphi }(\lambda ,\varphi ,\alpha ).`$
For oblate models with thin S-tube orbits ($`T_{\lambda \lambda }0`$, see §2.5.6), the analytical solution of (2.5.2) was derived by Bishop (1987) and by de Zeeuw & Hunter (1990). Robijn & de Zeeuw (1996) obtained the second-order velocity moments for models in which the thin tube orbits were thickened iteratively. Dejonghe & de Zeeuw (1988, Appendix D) found a general solution by integrating along characteristics. Evans (1990) gave an algorithm for solving (2.5.2) numerically, and Arnold (1995) computed a solution using characteristics without assuming a separable potential.

#### 2.5.3 Confocal elliptic coordinates: non-circular discs

In the principal plane $`z=0`$, the ellipsoidal coordinates reduce to confocal elliptic coordinates ($`\lambda ,\mu `$), with coordinate curves that are ellipses ($`\lambda `$) and hyperbolae ($`\mu `$), that share their foci on the symmetry $`y`$-axis. The potential of the perfect elliptic disc, with its surface density distribution stratified on concentric ellipses in the plane $`z=0`$ ($`\nu =\gamma `$), is of Stäckel form both in and outside this plane. By a superposition of perfect elliptic discs, one can construct other surface densities and corresponding disc potentials that are of Stäckel form in the plane $`z=0`$, but not necessarily outside it (Evans & de Zeeuw 1992). The expression for the potential in the disc is of the form (2.18) with $`g(\chi )0`$:

$$V_S=\frac{f(\lambda )f(\mu )}{\lambda \mu },$$
(2.24)
where again $`f(\tau )=(\tau +\alpha )G(\tau )`$, so that $`G(\tau )`$ equals the potential along the $`y`$-axis.

Omitting all terms with $`\nu `$ in (2.16), we obtain the Jeans equations for non-circular Stäckel discs

$`{\displaystyle \frac{T_{\lambda \lambda }}{\lambda }}+{\displaystyle \frac{T_{\lambda \lambda }T_{\mu \mu }}{2(\lambda \mu )}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\lambda }},`$
$`{\displaystyle \frac{T_{\mu \mu }}{\mu }}+{\displaystyle \frac{T_{\mu \mu }T_{\lambda \lambda }}{2(\mu \lambda )}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\mu }},`$
where now $`\rho `$ denotes a surface density. The parts of the $`y`$-axis between and beyond the foci are labeled by $`\lambda =\alpha `$ and $`\mu =\alpha `$, resulting in the continuity condition

$$T_{\lambda \lambda }(\alpha ,\alpha )=T_{\mu \mu }(\alpha ,\alpha ).$$
(2.26)
#### 2.5.4 Conical coordinates: scale-free triaxial limit

At large radii, the confocal ellipsoidal coordinates ($`\lambda ,\mu ,\nu `$) reduce to conical coordinates ($`r,\mu ,\nu `$), with $`r`$ the usual distance to the origin, i.e., $`r^2=x^2+y^2+z^2`$ and $`\mu `$ and $`\nu `$ angular coordinates on the sphere (Fig. 3). The potential $`V_S(r,\mu ,\nu )`$ is scale-free, and of the form

$$V_S=\stackrel{~}{F}(r)+\frac{F(\mu )F(\nu )}{r^2(\mu \nu )},$$
(2.27)
where $`\stackrel{~}{F}(r)`$ is arbitrary, and $`F(\tau )=(\tau +\alpha )(\tau +\gamma )G(\tau )`$, as in eq. (2.4).

The Jeans equations in conical coordinates follow from the general triaxial case (2.16) by going to large radii. Taking $`\lambda r^2\alpha \mu ,\nu `$, the stress components approach each other and we have

$$\frac{T_{\lambda \lambda }T_{\mu \mu }}{2(\lambda \mu )},\frac{T_{\lambda \lambda }T_{\nu \nu }}{2(\lambda \nu )}\frac{1}{r}0,\frac{}{\lambda }\frac{1}{2r}\frac{}{\lambda }.$$
(2.28)
Hence, after multiplying (2.16a) by $`2r`$, the Jeans equations for scale-free Stäckel models are

$`{\displaystyle \frac{T_{rr}}{r}}+{\displaystyle \frac{2T_{rr}T_{\mu \mu }T_{\nu \nu }}{r}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{r}},`$
$`{\displaystyle \frac{T_{\mu \mu }}{\mu }}+{\displaystyle \frac{T_{\mu \mu }T_{\nu \nu }}{2(\mu \nu )}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\mu }},`$ (2.29)
$`{\displaystyle \frac{T_{\nu \nu }}{\nu }}+{\displaystyle \frac{T_{\nu \nu }T_{\mu \mu }}{2(\nu \mu )}}`$ $`=`$ $`\rho {\displaystyle \frac{V_S}{\nu }}.`$
The general Jeans equations in conical coordinates, as derived by Evans et al. (2000), reduce to (2.5.4) for vanishing mixed second moments. At the transition points between the curves of constant $`\mu `$ and $`\nu `$ ($`\mu =\nu =\beta `$), the tensor components $`T_{\mu \mu }`$ and $`T_{\nu \nu }`$ coincide, resulting in the continuity condition

$$T_{\lambda \lambda }(r,\beta ,\beta )=T_{\varphi \varphi }(r,\beta ,\beta ).$$
(2.30)
#### 2.5.5 One-dimensional limits

There are several additional limiting cases with more symmetry for which the form of $`V_S`$ (Lynden–Bell 1962b) and the associated Jeans equations follow in a straightforward way from the expressions that were given above. We only mention spheres and circular discs.

When $`\alpha =\beta =\gamma `$, the variables $`\mu `$ and $`\nu `$ loose their meaning and the ellipsoidal coordinates reduce to spherical coordinates $`(r,\theta ,\varphi )`$. A steady-state spherical model without a preferred axis is invariant under a rotation over the angles $`\theta `$ and $`\varphi `$, so that we are left with only one Jeans equation in $`r`$, and $`T_{\theta \theta }=T_{\varphi \varphi }`$. This equation can readily be obtained from the CBE in spherical coordinates (e.g., Binney & Tremaine 1987). It also follows as a limit from the Jeans equations (2.16) for triaxial Stäckel models or from any of the above two-dimensional limiting cases. Consider for example the Jeans equations in conical coordinates (2.5.4), and take $`\mu \theta `$ and $`\nu \varphi `$. The stress components $`T_{rr}`$ and $`T_{\mu \mu }=T_{\nu \nu }=T_{\varphi \varphi }=T_{\theta \theta }`$ depend only $`r`$, so that we are left with

$$\frac{dT_{rr}}{dr}+\frac{2(T_{rr}T_{\theta \theta })}{r}=\rho \frac{dV_S}{dr},$$
(2.31)
which is the well-known result for non-rotating spherical systems (Binney & Tremaine 1987).

In a similar way, the one Jeans equation for the circular disc-case follows from, e.g., the first equation of (2.5.3) by taking $`\mu =\alpha `$ and replacing $`T_{\mu \mu }`$ by $`T_{\varphi \varphi }`$, where $`\varphi `$ is the azimuthal angle defined in §2.5.2. With $`\lambda +\alpha =R^2`$ this gives

$$\frac{dT_{RR}}{dR}+\frac{T_{RR}T_{\varphi \varphi }}{R}=\rho \frac{dV_S}{dR},$$
(2.32)
which may be compared with Binney & Tremaine (1987), their eq. (4.29).

#### 2.5.6 Thin tube orbits

Each of the three tube orbit families in a triaxial Stäckel model consists of a rotation in one of the ellipsoidal coordinates and oscillations in the other two (§2.2). The I-tubes, for example, rotate in $`\nu `$ and oscillate in $`\lambda `$ and $`\mu `$, with turning points $`\mu _1`$, $`\mu _2`$ and $`\lambda _0`$, so that a typical orbit fills the volume

$$\gamma \nu \beta ,\mu _1\mu \mu _2,\alpha \lambda \lambda _0.$$
(2.33)
When we restrict ourselves to infinitesimally thin I-tubes, i.e., $`\mu _1=\mu _2`$, there is no motion in the $`\mu `$-coordinate. Therefore, the second-order velocity moment in this coordinate is zero, and thus also the corresponding stress component $`T_{\mu \mu }^\mathrm{I}0`$. As a result, eq. (2.16b) reduces to an algebraic relation between $`T_{\lambda \lambda }^\mathrm{I}`$ and $`T_{\nu \nu }^\mathrm{I}`$. This relation can be used to eliminate $`T_{\nu \nu }^\mathrm{I}`$ and $`T_{\lambda \lambda }^\mathrm{I}`$ from the remaining Jeans equations (2.16a) and (2.16c) respectively.

HZ92 solved the resulting two first-order PDEs (their Appendix B) and showed that the same result is obtained by direct evaluation of the second-order velocity moments, using the thin I-tube DF. They derived similar solutions for thin O- and S-tubes, for which there is no motion in the $`\lambda `$-coordinate, so that $`T_{\lambda \lambda }^\mathrm{O}0`$ and $`T_{\lambda \lambda }^\mathrm{S}0`$, respectively.

In Stäckel discs we have – besides the flat box orbits – only one family of (flat) tube orbits. For infinitesimally thin tube orbits $`T_{\lambda \lambda }0`$, so that the Jeans equations (2.5.3) reduce to two different relations between $`T_{\mu \mu }`$ and the density and potential. In §3.4.4, we show how this places restrictions on the form of the density and we give the solution for $`T_{\mu \mu }`$. We also show that the general solution of (2.5.3), which we obtain in §3, contains the thin tube result. The same is true for the triaxial case: the general solution of (2.16), which we derive in §4, contains the three thin tube orbit solutions as special cases (§4.6.6).

### 2.6 All two-dimensional cases are similar

EL89 showed that the Jeans equations in oblate and prolate spheroidal coordinates, (2.5.1) and (2.5.2), can be transformed to a system that is equivalent to the two Jeans equations (2.5.3) in confocal elliptic coordinates. Evans et al. (2000) arrived at the same two-dimensional form for Stäckel models with a scale-free DF. We introduce a transformation which differs slightly from that of EL89, but has the advantage that it removes the singular denominators in the Jeans equations.

The Jeans equations (2.5.1) for prolate potentials can be simplified by introducing as dependent variables

$$𝒯_{\tau \tau }(\lambda ,\mu )=(\lambda +\beta )^{\frac{1}{2}}(\mu +\beta )^{\frac{1}{2}}(T_{\tau \tau }T_{\chi \chi }),\tau =\lambda ,\mu ,$$
(2.34)
so that the first two equations in (2.5.1) transform to

$`{\displaystyle \frac{𝒯_{\lambda \lambda }}{\lambda }}+{\displaystyle \frac{𝒯_{\lambda \lambda }𝒯_{\mu \mu }}{2(\lambda \mu )}}`$ $`=`$ $`(\lambda +\beta )^{\frac{1}{2}}(\mu +\beta )^{\frac{1}{2}}\left[\rho {\displaystyle \frac{V_S}{\lambda }}+{\displaystyle \frac{T_{\chi \chi }}{\lambda }}\right],`$
$`{\displaystyle \frac{𝒯_{\mu \mu }}{\mu }}+{\displaystyle \frac{𝒯_{\mu \mu }𝒯_{\lambda \lambda }}{2(\mu \lambda )}}`$ $`=`$ $`(\mu +\beta )^{\frac{1}{2}}(\lambda +\beta )^{\frac{1}{2}}\left[\rho {\displaystyle \frac{V_S}{\mu }}+{\displaystyle \frac{T_{\chi \chi }}{\mu }}\right].`$
The third Jeans equation (2.5.1) can be integrated in a straightforward fashion to give the $`\chi `$-dependence of $`T_{\chi \chi }`$. It is trivially satisfied for prolate models with $`g(\chi )0`$. Hence if, following EL89, we regard $`T_{\chi \chi }(\lambda ,\mu )`$ as a function which can be prescribed, then equations (2.6) have known right hand sides, and are therefore of the same form as those of the disc case (2.5.3). The singular denominator $`(\mu +\beta )`$ of (2.5.1) has disappeared, and there is a boundary condition

$$𝒯_{\mu \mu }(\lambda ,\beta )=0,$$
(2.36)
due to the second continuity condition of (2.5.1) and the definition (2.34).

A similar reduction applies for oblate potentials. The middle equation of (2.5.2) can be integrated to give the $`\varphi `$-dependence of $`T_{\varphi \varphi }`$, and is trivially satisfied for oblate models. The remaining two equations (2.5.2) transform to

$`{\displaystyle \frac{𝒯_{\lambda \lambda }}{\lambda }}+{\displaystyle \frac{𝒯_{\lambda \lambda }𝒯_{\nu \nu }}{2(\lambda \nu )}}`$ $`=`$ $`(\lambda +\alpha )^{\frac{1}{2}}(\alpha \nu )^{\frac{1}{2}}\left[\rho {\displaystyle \frac{V_S}{\lambda }}+{\displaystyle \frac{T_{\varphi \varphi }}{\lambda }}\right],`$
$`{\displaystyle \frac{𝒯_{\nu \nu }}{\nu }}+{\displaystyle \frac{𝒯_{\nu \nu }𝒯_{\lambda \lambda }}{2(\nu \lambda )}}`$ $`=`$ $`(\alpha \nu )^{\frac{1}{2}}(\lambda +\alpha )^{\frac{1}{2}}\left[\rho {\displaystyle \frac{V_S}{\nu }}+{\displaystyle \frac{T_{\varphi \varphi }}{\nu }}\right],`$
in terms of the dependent variables

$$𝒯_{\tau \tau }(\lambda ,\nu )=(\lambda +\alpha )^{\frac{1}{2}}(\alpha \nu )^{\frac{1}{2}}(T_{\tau \tau }T_{\varphi \varphi }),\tau =\lambda ,\nu .$$
(2.38)
We now have two boundary conditions

$$𝒯_{\lambda \lambda }(\alpha ,\nu )=0,𝒯_{\nu \nu }(\lambda ,\alpha )=0,$$
(2.39)
as a consequence of the last two continuity conditions of (2.5.2) and the definitions (2.38).

In the case of a scale-free DF, the stress components in the Jeans equations in conical coordinates (2.5.4) have the form $`T_{\tau \tau }=r^\zeta 𝒯_{\tau \tau }(\mu ,\nu )`$, with $`\zeta >0`$ and $`\tau =r,\mu ,\nu `$. After substitution and multiplication by $`r^{\zeta +1}`$, the first equation of (2.5.4) reduces to

$$(2\zeta )𝒯_{rr}+𝒯_{\mu \mu }+𝒯_{\nu \nu }=r^{\zeta +1}\rho \frac{V_S}{r}.$$
(2.40)
When $`\zeta =2`$, $`𝒯_{rr}`$ drops out, so that the relation between $`𝒯_{\mu \mu }`$ and $`𝒯_{\nu \nu }`$ is known and the remaining two Jeans equations can be readily solved (Evans et al. 2000). In all other cases, $`𝒯_{rr}`$ can be obtained from (2.40) once we have solved the last two equations of (2.5.4) for $`𝒯_{\mu \mu }`$ and $`𝒯_{\nu \nu }`$. This pair of equations is identical to the system of Jeans equations (2.5.3) for the case of disc potentials. The latter is the simplest form of the equivalent two-dimensional problem for all Stäckel models with two degrees of freedom. We solve it in the next section.

Once we have derived the solution of (2.5.3), we may obtain the solution for prolate Stäckel potentials by replacing all terms $`\rho V_s/\tau `$ $`(\tau =\lambda ,\mu )`$ by the right-hand side of (2.6) and substituting the transformations (2.34) for $`T_{\lambda \lambda }`$ and $`T_{\mu \mu }`$. Similarly, our unified notation makes the application of the solution of (2.5.3) to the oblate case and to models with a scale-free DF straightforward (§3.4).

## 3 The two-dimensional case

We first apply Riemann’s method to solve the Jeans equations (2.5.3) in confocal elliptic coordinates for Stäckel discs (§2.5.3). This involves finding a Riemann–Green function that describes the solution for a source point of stress. The full solution is then obtained in compact form by representing the known right-hand side terms as a sum of sources. In §3.2, we introduce an alternative approach, the singular solution method. Unlike Riemann’s method, this can be extended to the three-dimensional case, as we show in §4. We analyse the choice of the boundary conditions in detail in §3.3. In §3.4, we apply the two-dimensional solution to the axisymmetric and scale-free limits, and we also consider a Stäckel disc built with thin tube orbits.

### 3.1 Riemann’s method

After differentiating the first Jeans equation of (2.5.3) with respect to $`\mu `$ and eliminating terms in $`T_{\mu \mu }`$ by applying the second equation, we obtain a second-order partial differential equation (PDE) for $`T_{\lambda \lambda }`$ of the form

$$\frac{^2T_{\lambda \lambda }}{\lambda \mu }\frac{3}{2(\lambda \mu )}\frac{T_{\lambda \lambda }}{\lambda }+\frac{1}{2(\lambda \mu )}\frac{T_{\lambda \lambda }}{\mu }=U_{\lambda \lambda }(\lambda ,\mu ).$$
(3.1)
Here $`U_{\lambda \lambda }`$ is a known function given by

$$U_{\lambda \lambda }=\frac{1}{(\lambda \mu )^{\frac{3}{2}}}\frac{}{\mu }\left[(\lambda \mu )^{\frac{3}{2}}\rho \frac{V_S}{\lambda }\right]\frac{\rho }{2(\lambda \mu )}\frac{V_S}{\mu }.$$
(3.2)
We obtain a similar second-order PDE for $`T_{\mu \mu }`$ by interchanging $`\lambda \mu `$. Both PDEs can be solved by Riemann’s method. To solve them simultaneously, we define the linear second-order differential operator

$$=\frac{^2}{\lambda \mu }\frac{c_1}{\lambda \mu }\frac{}{\lambda }+\frac{c_2}{\lambda \mu }\frac{}{\mu },$$
(3.3)
with $`c_1`$ and $`c_2`$ constants to be specified. Hence, the more general second-order PDE

$$T=U,$$
(3.4)
with $`T`$ and $`U`$ functions of $`\lambda `$ and $`\mu `$ alone, reduces to those for the two stress components by taking

$`T=T_{\lambda \lambda }`$ $`:`$ $`c_1=\frac{3}{2},c_2=\frac{1}{2},U=U_{\lambda \lambda },`$
$`T=T_{\mu \mu }`$ $`:`$ $`c_1=\frac{1}{2},c_2=\frac{3}{2},U=U_{\mu \mu }.`$
In what follows, we introduce a Riemann–Green function $`𝒢`$ and incorporate the left-hand side of (3.4) into a divergence. Green’s theorem then allows us to rewrite the surface integral as a line integral over its closed boundary, which can be evaluated if $`𝒢`$ is chosen suitably. We determine the Riemann–Green function $`𝒢`$ which satisfies the required conditions, and then construct the solution.

#### 3.1.1 Application of Riemann’s method

We form a divergence by defining a linear operator $`^{}`$, called the adjoint of $``$ (e.g., Copson 1975), as

$$^{}=\frac{^2}{\lambda \mu }+\frac{}{\lambda }(\frac{c_1}{\lambda \mu })\frac{}{\mu }(\frac{c_2}{\lambda \mu }).$$
(3.6)
The combination $`𝒢TT^{}𝒢`$ is a divergence for any twice differentiable function $`𝒢`$ because

$$𝒢TT^{}𝒢=L/\lambda +M/\mu ,$$
(3.7)
where

$`L(\lambda ,\mu )`$ $`=`$ $`{\displaystyle \frac{𝒢}{2}}{\displaystyle \frac{T}{\mu }}{\displaystyle \frac{T}{2}}{\displaystyle \frac{𝒢}{\mu }}{\displaystyle \frac{c_1𝒢T}{\lambda \mu }},`$
$`M(\lambda ,\mu )`$ $`=`$ $`{\displaystyle \frac{𝒢}{2}}{\displaystyle \frac{T}{\lambda }}{\displaystyle \frac{T}{2}}{\displaystyle \frac{𝒢}{\lambda }}+{\displaystyle \frac{c_2𝒢T}{\lambda \mu }}.`$
We now apply the PDE (3.4) and the definition (3.6) in zero-subscripted variables $`\lambda _0`$ and $`\mu _0`$. We integrate the divergence (3.7) over the domain $`D=\{(\lambda _0,\mu _0)`$: $`\lambda \lambda _0\mathrm{},\mu \mu _0\alpha \}`$, with closed boundary $`\mathrm{\Gamma }`$ (Fig. 4). It follows by Green’s theorem that

$$\begin{array}{c}\underset{D}{}d\lambda _0d\mu _0\left(𝒢_0TT_0^{}𝒢\right)=\hfill \\ \hfill \underset{\mathrm{\Gamma }}{}d\mu _0L(\lambda _0,\mu _0)\underset{\mathrm{\Gamma }}{}d\lambda _0M(\lambda _0,\mu _0),\end{array}$$
(3.9)
where $`\mathrm{\Gamma }`$ is circumnavigated counter-clockwise. Here $`_0`$ and $`_0^{}`$ denote the operators (3.3) and (3.6) in zero-subscripted variables. We shall seek a Riemann–Green function $`𝒢(\lambda _0,\mu _0)`$ which solves the PDE

$$_0^{}𝒢=0,$$
(3.10)
in the interior of $`D`$. Then the left-hand side of (3.9) becomes $`_Dd\lambda _0d\mu _0𝒢(\lambda _0,\mu _0)U(\lambda _0,\mu _0)`$. The right-hand side of (3.9) has a contribution from each of the four sides of the rectangular boundary $`\mathrm{\Gamma }`$. We suppose that $`M(\lambda _0,\mu _0)`$ and $`L(\lambda _0,\mu _0)`$ decay sufficiently rapidly as $`\lambda _0\mathrm{}`$ so that the contribution from the boundary at $`\lambda _0=\mathrm{}`$ vanishes and the infinite integration over $`\lambda _0`$ converges. Partial integration of the remaining terms then gives for the boundary integral

$$\begin{array}{c}\underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0[(\frac{𝒢}{\lambda _0}\frac{c_2𝒢}{\lambda _0\mu _0})T\underset{\mu _0=\mu }{]}+\underset{\mu }{\overset{\alpha }{}}\mathrm{d}\mu _0[(\frac{𝒢}{\mu _0}+\frac{c_1𝒢}{\lambda _0\mu _0})T\underset{\lambda _0=\lambda }{]}\hfill \\ \hfill +\underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0[(\frac{T}{\lambda _0}+\frac{c_2T}{\lambda _0\mu _0})𝒢\underset{\mu _0=\alpha }{]}+𝒢(\lambda ,\mu )T(\lambda ,\mu ).\end{array}$$
(3.11)
We now impose on $`𝒢`$ the additional conditions

$$𝒢(\lambda ,\mu )=1,$$
(3.12)
and

$`{\displaystyle \frac{𝒢}{\lambda _0}}{\displaystyle \frac{c_2𝒢}{\lambda _0\mu _0}}=0`$ $`\mathrm{on}`$ $`\mu _0=\mu ,`$
$`{\displaystyle \frac{𝒢}{\mu _0}}+{\displaystyle \frac{c_1𝒢}{\lambda _0\mu _0}}=0`$ $`\mathrm{on}`$ $`\lambda _0=\lambda .`$
Then eq. (3.9) gives the explicit solution

$$\begin{array}{c}T(\lambda ,\mu )=\underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\mu }{\overset{\alpha }{}}d\mu _0𝒢(\lambda _0,\mu _0)U(\lambda _0,\mu _0)\hfill \\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0[(\frac{T}{\lambda _0}+\frac{c_2T}{\lambda _0\mu _0})𝒢\underset{\mu _0=\alpha }{]},\end{array}$$
(3.14)
for the stress component, once we have found the Riemann–Green function $`𝒢`$.

#### 3.1.2 The Riemann–Green function

Our prescription for the Riemann–Green function $`𝒢(\lambda _0,\mu _0)`$ is that it satisfies the PDE (3.10) as a function of $`\lambda _0`$ and $`\mu _0`$, and that it satisfies the boundary conditions (3.12) and (3.1.1) at the specific values $`\lambda _0=\lambda `$ and $`\mu _0=\mu `$. Consequently $`𝒢`$ depends on two sets of coordinates. Henceforth, we denote it as $`𝒢(\lambda ,\mu ;\lambda _0,\mu _0)`$.

An explicit expression for the Riemann–Green function which solves (3.10) is (Copson 1975)

$$𝒢(\lambda ,\mu ;\lambda _0,\mu _0)=\frac{(\lambda _0\mu _0)^{c_2}(\lambda \mu _0)^{c_1c_2}}{(\lambda \mu )^{c_1}}F(w),$$
(3.15)
where the parameter $`w`$ is defined as

$$w=\frac{(\lambda _0\lambda )(\mu _0\mu )}{(\lambda _0\mu _0)(\lambda \mu )},$$
(3.16)
and $`F(w)`$ is to be determined. Since $`w=0`$ when $`\lambda _0=\lambda `$ or $`\mu _0=\mu `$, it follows from (3.12) that the function $`F`$ has to satisfy $`F(0)=1`$. It is straightforward to verify that $`𝒢`$ satisfies the conditions (3.1.1), and that eq. (3.10) reduces to the following ordinary differential equation for $`F(w)`$

$$w(1w)F^{\prime \prime }+[1(2+c_1c_2)w]F^{}c_1(1c_2)F=0.$$
(3.17)
This is a hypergeometric equation (e.g., Abramowitz & Stegun 1965), and its unique solution satisfying $`F(0)=1`$ is

$$F(w)={}_{2}{}^{}F_{1}^{}(c_1,1c_2;1;w).$$
(3.18)
The Riemann–Green function (3.15) represents the influence at a field point at $`(\lambda ,\mu )`$ due to a source point at $`(\lambda _0,\mu _0)`$. Hence it satisfies the PDE

$$𝒢(\lambda ,\mu ;\lambda _0,\mu _0)=\delta (\lambda _0\lambda )\delta (\mu _0\mu ).$$
(3.19)
The first right-hand side term of the solution (3.14) is a sum over the sources in $`D`$ which are due to the inhomogeneous term $`U`$ in the PDE (3.4). That PDE is hyperbolic with characteristic variables $`\lambda `$ and $`\mu `$. By choosing to apply Green’s theorem to the domain $`D`$, we made it the domain of dependence (Strauss 1992) of the field point $`(\lambda ,\mu )`$ for (3.4), and hence we implicitly decided to integrate that PDE in the direction of decreasing $`\lambda `$ and decreasing $`\mu `$.

The second right-hand side term of the solution (3.14) represents the solution to the homogeneous PDE $`T=0`$ due to the boundary values of $`T`$ on the part of the boundary $`\mu =\alpha `$ which lies within the domain of dependence. There is only one boundary term because we implicitly require that $`T(\lambda ,\mu )0`$ as $`\lambda \mathrm{}`$. We verify in §3.1.4 that this requirement is indeed satisfied.

#### 3.1.3 The disc solution

We obtain the Riemann–Green functions for $`T_{\lambda \lambda }`$ and $`T_{\mu \mu }`$, labeled as $`𝒢_{\lambda \lambda }`$ and $`𝒢_{\mu \mu }`$, respectively, from expressions (3.15) and (3.18) by substitution of the values for the constants $`c_1`$ and $`c_2`$ from (3.1). The hypergeometric function in $`𝒢_{\lambda \lambda }`$ is the complete elliptic integral of the second kind222We use the definition $`E(w)=_0^{\frac{\pi }{2}}d\theta \sqrt{1w\mathrm{sin}^2\theta }`$, $`E(w)`$. The hypergeometric function in $`𝒢_{\mu \mu }`$ can also be expressed in terms of $`E(w)`$ using eq. (15.2.15) of Abramowitz & Stegun (1965), so that we can write

$$𝒢_{\lambda \lambda }(\lambda ,\mu ;\lambda _0,\mu _0)=\frac{(\lambda _0\mu _0)^{\frac{3}{2}}}{(\lambda \mu )^{\frac{1}{2}}}\frac{2E(w)}{\pi (\lambda _0\mu )},$$
(3.20a)
$$𝒢_{\mu \mu }(\lambda ,\mu ;\lambda _0,\mu _0)=\frac{(\lambda _0\mu _0)^{\frac{3}{2}}}{(\lambda \mu )^{\frac{1}{2}}}\frac{2E(w)}{\pi (\lambda \mu _0)},$$
(3.20b)
Substituting these into (3.14) gives the solution of the stress components throughout the disc as

$$\begin{array}{c}T_{\lambda \lambda }(\lambda ,\mu )=\frac{2}{\pi (\lambda \mu )^{\frac{1}{2}}}\{\hfill \\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\mu }{\overset{\alpha }{}}d\mu _0\frac{E(w)}{(\lambda _0\mu )}\{\frac{}{\mu _0}\left[(\lambda _0\mu _0)^{\frac{3}{2}}\rho \frac{V_S}{\lambda _0}\right]\frac{(\lambda _0\mu _0)^{\frac{1}{2}}}{2}\rho \frac{V_S}{\mu _0}\}\\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0[\frac{E(w)}{(\lambda _0\mu )}\underset{\mu _0=\alpha }{]}(\lambda _0+\alpha )\frac{\mathrm{d}}{\mathrm{d}\lambda _0}[(\lambda _0+\alpha )^{\frac{1}{2}}T_{\lambda \lambda }(\lambda _0,\alpha )]\},\end{array}$$
(3.21a)
$$\begin{array}{c}T_{\mu \mu }(\lambda ,\mu )=\frac{2}{\pi (\lambda \mu )^{\frac{1}{2}}}\{\hfill \\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\mu }{\overset{\alpha }{}}d\mu _0\frac{E(w)}{(\lambda \mu _0)}\{\frac{}{\lambda _0}\left[(\lambda _0\mu _0)^{\frac{3}{2}}\rho \frac{V_S}{\mu _0}\right]+\frac{(\lambda _0\mu _0)^{\frac{1}{2}}}{2}\rho \frac{V_S}{\lambda _0}\}\\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0[\frac{E(w)}{(\lambda \mu _0)}\underset{\mu _0=\alpha }{]}\frac{\mathrm{d}}{\mathrm{d}\lambda _0}[(\lambda _0+\alpha )^{\frac{3}{2}}T_{\mu \mu }(\lambda _0,\alpha )]\}.\end{array}$$
(3.21b)
This solution depends on $`\rho `$ and $`V_S`$, which are assumed to be known, and on $`T_{\lambda \lambda }(\lambda ,\alpha )`$ and $`T_{\mu \mu }(\lambda ,\alpha )`$, i.e., the stress components on the part of the $`y`$-axis beyond the foci. Because these two stress components satisfy the first Jeans equation of (2.5.3) at $`\mu =\alpha `$, we are only free to choose one of them, say $`T_{\mu \mu }(\lambda ,\alpha )`$. $`T_{\lambda \lambda }(\lambda ,\alpha )`$ then follows by integrating this first Jeans equation with respect to $`\lambda `$, using the continuity condition (2.26) and requiring that $`T_{\lambda \lambda }(\lambda ,\alpha )0`$ as $`\lambda \mathrm{}`$.

#### 3.1.4 Consistency check

We now investigate the behaviour of our solutions at large distances and verify that our working hypothesis concerning the radial fall-off of the functions $`L`$ and $`M`$ in eq. (3.1.1) is correct. The solution (3.14) consists of two components: an area integral due to the inhomogeneous right-hand side term of the PDE (3.4), and a single integral due to the boundary values. We examine them in turn to obtain the conditions for the integrals to converge. Next, we parameterise the behaviour of the density and potential at large distances and apply it to the solution (3.21) and to the energy equation (2.10) to check if the convergence conditions are satisfied for physical potential-density pairs.

As $`\lambda _0\mathrm{}`$, $`w`$ tends to the finite limit $`(\mu _0\mu )/(\lambda \mu )`$. Hence $`E(w)`$ is finite, and so, by (3.20), $`𝒢_{\lambda \lambda }=𝒪(\lambda _0^{1/2})`$ and $`𝒢_{\mu \mu }=𝒪(\lambda _0^{3/2})`$. Suppose now that $`U_{\lambda \lambda }(\lambda _0,\mu _0)=𝒪(\lambda _0^{l_11})`$ and $`U_{\mu \mu }(\lambda _0,\mu _0)=𝒪(\lambda _0^{m_11})`$ as $`\lambda _0\mathrm{}`$. The area integrals in the solution (3.14) then converge, provided that $`l_1>\frac{1}{2}`$ and $`m_1>\frac{3}{2}`$. These requirements place restrictions on the behaviour of the density $`\rho `$ and potential $`V_S`$ which we examine below. Since $`𝒢_{\lambda \lambda }(\lambda ,\mu ;\lambda _0,\mu _0)`$ is $`𝒪(\lambda ^{1/2})`$ as $`\lambda \mathrm{}`$, the area integral component of $`T_{\lambda \lambda }(\lambda ,\mu )`$ behaves as $`𝒪(\lambda ^{1/2}_\lambda ^{\mathrm{}}\lambda _0^{l_11/2}𝑑\lambda _0)`$ and so is $`𝒪(\lambda ^{l_1})`$. Similarly, with $`𝒢_{\mu \mu }(\lambda ,\mu ;\lambda _0,\mu _0)=𝒪(\lambda ^{3/2})`$ as $`\lambda \mathrm{}`$, the first component of $`T_{\mu \mu }(\lambda ,\mu )`$ is $`𝒪(\lambda _0^{m_1})`$.

To analyse the second component of the solution (3.14), we suppose that the boundary value $`T_{\lambda \lambda }(\lambda _0,\alpha )=𝒪(\lambda _0^{l_2})`$ and $`T_{\mu \mu }(\lambda _0,\alpha )=𝒪(\lambda _0^{m_2})`$ as $`\lambda _0\mathrm{}`$. A similar analysis then shows that the boundary integrals converge, provided that $`l_2>\frac{1}{2}`$ and $`m_2>\frac{3}{2}`$, and that the second components of $`T_{\lambda \lambda }(\lambda ,\mu )`$ and $`T_{\mu \mu }(\lambda ,\mu )`$ are $`𝒪(\lambda ^{l_2})`$ and $`𝒪(\lambda ^{m_2})`$ as $`\lambda \mathrm{}`$, respectively.

We conclude that the convergence of the integrals in the solution (3.14) requires that $`T_{\lambda \lambda }(\lambda ,\mu )`$ and $`T_{\mu \mu }(\lambda ,\mu )`$ decay at large distance as $`𝒪(\lambda ^l)`$ with $`l>\frac{1}{2}`$ and $`𝒪(\lambda ^m)`$ with $`m>\frac{3}{2}`$, respectively. The requirements which we have imposed on $`U(\lambda _0,\mu _0)`$ and $`T(\lambda _0,\alpha )`$ cause the contributions to $`_\mathrm{\Gamma }𝑑\mu _0L(\lambda _0,\mu _0)`$ in Green’s formula (3.9) from the segment of the path at large $`\lambda _0`$ to be negligible in all cases.

Having obtained the requirements for the Riemann–Green function analysis to be valid, we now investigate the circumstances in which they apply. Following Arnold et al. (1994), we consider densities $`\rho `$ that decay as $`N(\mu )\lambda ^{s/2}`$ at large distances. We suppose that the function $`G(\tau )`$ introduced in eq. (2.4) is $`𝒪(\tau ^\delta )`$ for $`\frac{1}{2}\delta <0`$ as $`\tau \mathrm{}`$. The lower limit $`\delta =\frac{1}{2}`$ corresponds to a potential due to a finite total mass, while the upper limit restricts it to potentials that decay to zero at large distances.

For the disc potential (2.24), we then have that $`f(\tau )=𝒪(\tau ^{\delta +1})`$ when $`\tau \mathrm{}`$. Using the definition (3.2), we obtain

$$U_{\lambda \lambda }(\lambda ,\mu )=\frac{f^{}(\mu )f^{}(\lambda )}{2(\lambda \mu )^2}\rho +\frac{V_S+f^{}(\lambda )}{(\lambda \mu )}\frac{\rho }{\mu },$$
(3.22a)
$$U_{\mu \mu }(\lambda ,\mu )=\frac{f^{}(\lambda )f^{}(\mu )}{2(\lambda \mu )^2}\rho \frac{V_S+f^{}(\mu )}{(\lambda \mu )}\frac{\rho }{\lambda },$$
(3.22b)
where $`\rho `$ is the surface density of the disc. It follows that $`U_{\lambda \lambda }(\lambda ,\mu )`$ is generally the larger and is $`𝒪(\lambda ^{\delta s/21})`$ as $`\lambda \mathrm{}`$, whereas $`U_{\mu \mu }(\lambda ,\mu )`$ is $`𝒪(\lambda ^{2s/2})`$. Hence, for the components of the stresses (3.21) we have $`T_{\lambda \lambda }=𝒪(\lambda ^{\delta s/2})`$ and $`T_{\mu \mu }=𝒪(\lambda ^{1s/2})`$. This estimate for $`U_{\lambda \lambda }`$ assumes that $`\rho /\mu `$ is also $`𝒪(\lambda ^{s/2})`$. It is too high if the density becomes independent of angle at large distances, as it does for discs with $`s<3`$ (Evans & de Zeeuw 1992). Using these estimates with the requirements for integral convergence that were obtained earlier, we obtain the conditions $`s>2\delta +1`$ and $`s>1`$, respectively, for inhomogeneous terms in $`T_{\lambda \lambda }(\lambda ,\mu )`$ and $`T_{\mu \mu }(\lambda ,\mu )`$ to be valid solutions. The second condition implies the first because $`\delta <0`$.

With $`V_S(\lambda ,\mu )=𝒪(\lambda ^\delta )`$ at large $`\lambda `$, it follows from the energy equation (2.10) for bound orbits that the second-order velocity moments $`v_\tau ^2`$ cannot exceed $`𝒪(\lambda ^\delta )`$, and hence that stresses $`T_{\tau \tau }=\rho v_\tau ^2`$ cannot exceed $`𝒪(\lambda ^{\delta s/2})`$. This implies for $`T_{\lambda \lambda }(\lambda ,\mu )`$ that $`s>2\delta +1`$, and for $`T_{\mu \mu }(\lambda ,\mu )`$ we have the more stringent requirement that $`s>2\delta +3`$. This last requirement is unnecessarily restrictive, but an alternative form of the solution is needed to do better. Since that alternative form arises naturally with the singular solution method, we return to this issue in §3.2.6.

Thus, for the Riemann–Green solution to apply, we find the conditions $`s>1`$ and $`\frac{1}{2}\delta <0`$. These conditions are satisfied for the perfect elliptic disk $`(s=3,\delta =\frac{1}{2})`$, and for many other separable discs (Evans & de Zeeuw 1992).

#### 3.1.5 Relation to the EL89 analysis

EL89 solve for the difference $`\mathrm{\Delta }T_{\lambda \lambda }T_{\mu \mu }`$ using a Green’s function method which is essentially equivalent to the approach used here. EL89 give the Fourier transform of their Green’s function, but do not invert it. We give the Riemann–Green function for $`\mathrm{\Delta }`$ in Appendix A, and then rederive it by a Laplace transform analysis. Our Laplace transform analysis can be recast in terms of Fourier transforms. When we do this, we obtain a result which differs from that of EL89.

### 3.2 Singular Solution Superposition

We have solved the disc problem (2.5.3) by combining the two Jeans equations into a single second-order PDE in one of the stress components, and then applying Riemann’s method to it. However, Riemann’s method and other standard techniques do not carry over to a single third-order PDE in one dependent variable, which is the best that one could expect to have in the general case. We therefore introduce an alternative but equivalent method of solution, also based on the superposition of source points. In constrast to Riemann’s method, this singular solution method is applicable to the general case of triaxial Stäckel models.

#### 3.2.1 Simplified Jeans equations

We define new independent variables

$`S_{\lambda \lambda }(\lambda ,\mu )`$ $`=`$ $`|\lambda \mu |^{\frac{1}{2}}T_{\lambda \lambda }(\lambda ,\mu ),`$
$`S_{\mu \mu }(\lambda ,\mu )`$ $`=`$ $`|\mu \lambda |^{\frac{1}{2}}T_{\mu \mu }(\lambda ,\mu ),`$
where $`|.|`$ denotes absolute value, introduced to make the square root single-valued with respect to cyclic permutation of $`\lambda \mu \lambda `$. The Jeans equations (2.5.3) can then be written in the form

$$\frac{S_{\lambda \lambda }}{\lambda }\frac{S_{\mu \mu }}{2(\lambda \mu )}=|\lambda \mu |^{\frac{1}{2}}\rho \frac{V_S}{\lambda }g_1(\lambda ,\mu ),$$
(3.24a)
$$\frac{S_{\mu \mu }}{\mu }\frac{S_{\lambda \lambda }}{2(\mu \lambda )}=|\mu \lambda |^{\frac{1}{2}}\rho \frac{V_S}{\mu }g_2(\lambda ,\mu ).$$
(3.24b)
For given density and potential, $`g_1`$ and $`g_2`$ are known functions of $`\lambda `$ and $`\mu `$. Next, we consider a simplified form of (3.24) by taking for $`g_1`$ and $`g_2`$, respectively

$$\stackrel{~}{g}_1(\lambda ,\mu )=0,\stackrel{~}{g}_2(\lambda ,\mu )=\delta (\lambda _0\lambda )\delta (\mu _0\mu ),$$
(3.25)
with $`\beta \mu \mu _0\alpha \lambda \lambda _0`$. A similar set of simplified equations is obtained by interchanging the expressions for $`\stackrel{~}{g}_1`$ and $`\stackrel{~}{g}_2`$. We refer to solutions of these simplified Jeans equations as singular solutions.

Singular solutions can be interpreted as contributions to the stresses at a fixed point $`(\lambda ,\mu )`$ due to a source point in $`(\lambda _0,\mu _0)`$ (Fig. 4). The full stress at the field point can be obtained by adding all source point contributions, each with a weight that depends on the local density and potential. In what follows, we derive the singular solutions, and then use this superposition principle to construct the solution for the Stäckel discs in §3.2.6.

#### 3.2.2 Homogeneous boundary problem

The choice (3.25) places constraints on the functional form of $`S_{\lambda \lambda }`$ and $`S_{\mu \mu }`$. The presence of the delta-functions in $`\stackrel{~}{g}_2`$ requires that $`S_{\mu \mu }`$ contains a term $`\delta (\lambda _0\lambda )(\mu _0\mu )`$, with the step-function

$$(xx_0)=\{\begin{array}{cc}0,\hfill & x<x_0,\hfill \\ 1,\hfill & xx_0.\hfill \end{array}$$
(3.26)
Since $`^{}(y)=\delta (y)`$, it follows that, by taking the partial derivative of $`\delta (\lambda _0\lambda )(\mu _0\mu )`$ with respect to $`\mu `$, the delta-functions are balanced. There is no balance when $`S_{\lambda \lambda }`$ contains $`\delta (\lambda _0\lambda )`$, and similarly neither stress components can contain $`\delta (\mu _0\mu )`$. We can, however, add a function of $`\lambda `$ and $`\mu `$ to both components, multiplied by $`(\lambda _0\lambda )(\mu _0\mu )`$. In this way, we obtain a singular solution of the form

$`S_{\lambda \lambda }`$ $`=`$ $`A(\lambda ,\mu )(\lambda _0\lambda )(\mu _0\mu ),`$
$`S_{\mu \mu }`$ $`=`$ $`B(\lambda ,\mu )(\lambda _0\lambda )(\mu _0\mu )\delta (\lambda _0\lambda )(\mu _0\mu ),`$
in terms of functions $`A`$ and $`B`$ that have to be determined. Substituting these forms in the simplified Jeans equations and matching terms yields two homogeneous equations

$$\frac{A}{\lambda }\frac{B}{2(\lambda \mu )}=0,\frac{B}{\mu }\frac{A}{2(\mu \lambda )}=0,$$
(3.28)
and two boundary conditions

$$A(\lambda _0,\mu )=\frac{1}{2(\lambda _0\mu )},B(\lambda ,\mu _0)=0.$$
(3.29)
Two alternative boundary conditions which are useful below can be found as follows. Integrating the first of the equations (3.28) with respect to $`\lambda `$ on $`\mu =\mu _0`$, where $`B(\lambda ,\mu _0)=0`$, gives

$$A(\lambda ,\mu _0)=\frac{1}{2(\lambda _0\mu _0)}.$$
(3.30)
Similarly, integrating the second of equations (3.28) with respect to $`\mu `$ on $`\lambda =\lambda _0`$ where $`A`$ is known gives

$$B(\lambda _0,\mu )=\frac{\mu _0\mu }{4(\lambda _0\mu _0)(\lambda _0\mu )}.$$
(3.31)
Even though expressions (3.30) and (3.31) do not add new information, they will be useful for identifying contour integral formulas in the analysis which follows.

We have reduced the problem of solving the Jeans equations (2.5.3) for Stäckel discs to a two-dimensional boundary problem. We solve this problem by first deriving a one-parameter particular solution (§3.2.3) and then making a linear combination of particular solutions with different values of their free parameter, such that the four boundary expressions are satisfied simultaneously (§3.2.4). This gives the solution of the homogeneous boundary problem.

#### 3.2.3 Particular solution

To find a particular solution of the homogeneous equations (3.28) with one free parameter $`z`$, we take as an Ansatz

$`A(\lambda ,\mu )`$ $``$ $`(\lambda \mu )^{a_1}(z\lambda )^{a_2}(z\mu )^{a_3},`$
$`B(\lambda ,\mu )`$ $``$ $`(\lambda \mu )^{b_1}(z\lambda )^{b_2}(z\mu )^{b_3},`$
with $`a_i`$ and $`b_i`$ $`(i=1,2,3)`$ all constants. Hence,

$`{\displaystyle \frac{A}{\lambda }}`$ $`=`$ $`A\left({\displaystyle \frac{a_1}{\lambda \mu }}{\displaystyle \frac{a_2}{z\lambda }}\right)={\displaystyle \frac{1}{2(\lambda \mu )}}\left(2a_1A{\displaystyle \frac{z\mu }{z\lambda }}\right),`$
$`{\displaystyle \frac{B}{\mu }}`$ $`=`$ $`B\left({\displaystyle \frac{b_1}{\mu \lambda }}{\displaystyle \frac{b_3}{z\mu }}\right)={\displaystyle \frac{1}{2(\mu \lambda )}}\left(2b_1B{\displaystyle \frac{z\lambda }{z\mu }}\right),`$
where we have set $`a_2=a_1`$ and $`b_3=b_1`$. Taking $`a_1=b_1=\frac{1}{2}`$, the homogeneous equations are satisfied if

$$\frac{z\lambda }{z\mu }=\frac{A}{B}=\frac{(z\lambda )^{\frac{1}{2}b_2}}{(z\mu )^{\frac{1}{2}a_3}},$$
(3.34)
so, $`a_3=b_2=\frac{3}{2}`$. We denote the resulting solutions as

$$A^P(\lambda ,\mu )=\frac{|\lambda \mu |^{\frac{1}{2}}}{(z\lambda )^{\frac{1}{2}}(z\mu )^{\frac{3}{2}}},$$
(3.35a)
$$B^P(\lambda ,\mu )=\frac{|\mu \lambda |^{\frac{1}{2}}}{(z\mu )^{\frac{1}{2}}(z\lambda )^{\frac{3}{2}}}.$$
(3.35b)
These particular solutions follow from each other by cyclic permutation $`\lambda \mu \lambda `$, as is required from the symmetry of the homogeneous equations (3.28).

#### 3.2.4 The homogeneous solution

We now consider a linear combination of the particular solution (3.35) by integrating it over the free parameter $`z`$, which we assume to be complex. We choose the integration contours in the complex $`z`$-plane, such that the four boundary expressions can be satisfied simultaneously.

We multiply $`B^P(\lambda ,\mu )`$ by $`(z\mu _0)^{\frac{1}{2}}`$, and integrate it over the closed contour $`C^\mu `$ (Fig. 5). When $`\mu =\mu _0`$, the integrand is analytic within $`C^\mu `$, so that the integral vanishes by Cauchy’s theorem. Since both the multiplication factor and the integration are independent of $`\lambda `$ and $`\mu `$, it follows from the superposition principle that the homogeneous equations are still satisfied. In this way, the second of the boundary expressions (3.29) is satisfied.

Next, we also multiply $`B^P(\lambda ,\mu )`$ by $`(z\lambda _0)^{\frac{1}{2}}`$, so that the contour $`C^\lambda `$ (Fig. 5) encloses a double pole when $`\lambda =\lambda _0`$. From the Residue theorem (e.g., Conway 1973), it then follows that

$$\begin{array}{c}\underset{C^\lambda }{}\frac{(z\mu _0)^{\frac{1}{2}}}{(z\lambda _0)^{\frac{1}{2}}}B^P(\lambda _0,\mu )dz=\underset{C^\lambda }{}\frac{(z\mu _0)^{\frac{1}{2}}(\lambda _0\mu )^{\frac{1}{2}}}{(z\mu )^{\frac{1}{2}}(z\lambda _0)^2}dz=\hfill \\ \hfill 2\pi i(\lambda _0\mu )^{\frac{1}{2}}[\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{z\mu _0}{z\mu }\right)^{\frac{1}{2}}\underset{z=\lambda _0}{]}=\frac{\pi i(\mu _0\mu )}{(\lambda _0\mu _0)^{\frac{1}{2}}(\lambda _0\mu )},\end{array}$$
(3.36)
which equals the boundary expression (3.31), up to the factor $`4\pi i(\lambda _0\mu _0)^{\frac{1}{2}}`$.

Taking into account the latter factor, and the ratio (3.34) of $`A`$ and $`B`$, we postulate as homogeneous solution

$$A(\lambda ,\mu )=\frac{1}{4\pi i}\frac{|\lambda \mu |^{\frac{1}{2}}}{|\lambda _0\mu _0|^{\frac{1}{2}}}\underset{C}{}\frac{(z\mu _0)^{\frac{1}{2}}\mathrm{d}z}{(z\lambda )^{\frac{1}{2}}(z\mu )^{\frac{3}{2}}(z\lambda _0)^{\frac{1}{2}}},$$
(3.37a)
$$B(\lambda ,\mu )=\frac{1}{4\pi i}\frac{|\mu \lambda |^{\frac{1}{2}}}{|\lambda _0\mu _0|^{\frac{1}{2}}}\underset{C}{}\frac{(z\mu _0)^{\frac{1}{2}}\mathrm{d}z}{(z\mu )^{\frac{1}{2}}(z\lambda )^{\frac{3}{2}}(z\lambda _0)^{\frac{1}{2}}},$$
(3.37b)
with the choice for the contour $`C`$ still to be specified.

The integrands in (3.37) consist of multi-valued functions that all come in pairs $`(z\tau )^{1/2m}(z\tau _0)^{1/2n}`$, for integer $`m`$ and $`n`$, and for $`\tau `$ being either $`\lambda `$ or $`\mu `$. Hence, we can make the integrands single-valued by specifying two cuts in the complex $`z`$-plane, one from $`\mu `$ to $`\mu _0`$ and one from $`\lambda `$ to $`\lambda _0`$. The integrands are now analytic in the cut plane away from its cuts and behave as $`z^2`$ at large distances, so that the integral over a circular contour with infinite radius is zero333We evaluate the square roots as $`(z\tau )^{\frac{1}{2}}=|z\tau |\mathrm{exp}i\mathrm{arg}(z\tau )`$ with $`|\mathrm{arg}(z\tau )|\pi `$.. Connecting the simple contours $`C^\lambda `$ and $`C^\mu `$ with this circular contour shows that the cumulative contribution from each of these contours cancels. As a consequence, every time we integrate over the contour $`C^\lambda `$, we will obtain the same result by integrating over $`C^\mu `$ instead. This means we integrate over $`C^\mu `$ and take the negative of the result or, equally, integrate over $`C^\mu `$ in clockwise direction.

For example, we obtained the boundary expression for $`B`$ in (3.36) by applying the Residue theorem to the double pole enclosed by the contour $`C^\lambda `$. The evaluation of the integral becomes less straightforward when we consider the contour $`C^\mu `$ instead. Wrapping the contour around the branch points $`\mu `$ and $`\mu _0`$ (Fig. 6), one may easily verify that the contribution from the two arcs vanishes if their radius goes to zero. Taking into account the change in phase when going around the two branch points, one may show that the contributions from the two remaining parts of the contour, parallel to the real axis, are equivalent. Hence, we arrive at the following (real) integral

$$B(\lambda _0,\mu )=\frac{1}{2\pi }\frac{(\lambda \mu _0)^{\frac{1}{2}}}{(\lambda _0\mu _0)^{\frac{1}{2}}}\underset{\mu }{\overset{\mu _0}{}}\frac{\mathrm{d}t}{(\lambda _0t)^2}\sqrt{\frac{\mu _0t}{t\mu }}.$$
(3.38)
The substitution

$$t=\mu +\frac{(\mu _0\mu )(\lambda _0\mu _0)\mathrm{sin}^2\theta }{(\mu _0\mu )s(\lambda _0\mu )}$$
(3.39)
then indeed gives the correct boundary expression (3.31).

When we take $`\mu =\mu _0`$ in (3.37b), we are left with the integrand $`(z\lambda )^{3/2}(z\lambda _0)^{1/2}`$. This is analytic within the contour $`C^\mu `$ and hence it follows from Cauchy’s theorem that there is no contribution. However, if we take the contour $`C^\lambda `$ instead, it is not clear at once that the integral indeed is zero. To evaluate the complex integral we wrap the contour $`C^\lambda `$ around the branch points $`\lambda `$ and $`\lambda _0`$ (Fig. 6). There will be no contribution from the arc around $`\lambda _0`$ if its radius goes to zero. However, since the integrand involves the term $`z\lambda `$ with power $`\frac{3}{2}`$, the contribution from the arc around $`\lambda `$ is of the order $`ϵ^{1/2}`$ and hence goes to infinity if its radius $`ϵ>0`$ reduces to zero. If we let the two remaining straight parts of the contour run from $`\lambda +ϵ`$ to $`\lambda _0`$, then their cumulative contribution becomes proportional to $`\mathrm{tan}\theta (ϵ)`$, with $`\theta (ϵ)`$ approaching $`\frac{\pi }{2}`$ when $`ϵ`$ reduces to zero. Hence, both the latter contribution and the contribution from the arc around $`\lambda `$ approaches infinity. However, careful investigation of their limiting behaviour shows that they cancel when $`ϵ`$ reaches zero, as is required for the boundary expression $`B(\lambda ,\mu _0)=0`$.

We have shown that the use of $`C^\lambda `$ and $`C^\mu `$ gives the same result, but the effort to evaluate the contour integral varies between the two choices. The boundary expressions for $`A(\lambda ,\mu )`$, (3.29) and (3.30) are obtained most easily if we consider $`C^\lambda `$ when $`\lambda =\lambda _0`$ and $`C^\mu `$ when $`\mu =\mu _0`$. In both cases the integrand in (3.37a) has a single pole within the chosen contour, so that the boundary expressions follow by straightforward application of the Residue theorem.

We now have proven that the homogeneous solution (3.37) solves the homogeneous equations (3.28), satisfies the boundary values (3.29)–(3.31) separately and, from the observation that $`C^\lambda `$ and $`C^\mu `$ produce the same result, also simultaneously.

#### 3.2.5 Evaluation of the homogeneous solution

The homogeneous solution (3.37) consists of complex contour integrals, which we transform to real integrals by wrapping the contours $`C^\lambda `$ and $`C^\mu `$ around the corresponding pair of branch points (Fig. 6). To have no contribution from the arcs around the branch points, we choose the (combination of) contours such that the terms in the integrand involving these branch points have powers larger than $`1`$. In this way, we can always evaluate the complex integral as a (real) integral running from one branch point to the other.

In the homogeneous solution (3.37a) for $`A`$ we choose $`C=C^\lambda `$ and in (3.37b) for $`B`$ we take $`C=C^\mu `$. Taking into account the changes in phase when going around the branch points, we obtain the following expressions for the homogeneous solution

$$A(\lambda ,\mu )=\frac{1}{2\pi }\frac{|\lambda \mu |^{\frac{1}{2}}}{|\lambda _0\mu _0|^{\frac{1}{2}}}\underset{\lambda }{\overset{\lambda _0}{}}\frac{\mathrm{d}t}{t\mu }\sqrt{\frac{t\mu _0}{(t\lambda )(t\mu )(\lambda _0t)}},$$
(3.40a)
$$B(\lambda ,\mu )=\frac{1}{2\pi }\frac{|\lambda \mu |^{\frac{1}{2}}}{|\lambda _0\mu _0|^{\frac{1}{2}}}\underset{\mu }{\overset{\mu _0}{}}\frac{\mathrm{d}t}{\lambda t}\sqrt{\frac{\mu _0t}{(\lambda t)(t\mu )(\lambda _0t)}}.$$
(3.40b)
By a parameterisation of the form (3.39), or by using an integral table (e.g., Byrd & Friedman 1971), expressions (3.40) can be written conveniently in terms of the complete elliptic integral of the second kind, $`E`$, and its derivative $`E^{}`$

$$A(\lambda ,\mu ;\lambda _0,\mu _0)=\frac{E(w)}{\pi (\lambda _0\mu )},$$
(3.41a)
$$B(\lambda ,\mu ;\lambda _0,\mu _0)=\frac{2wE^{}(w)}{\pi (\lambda _0\lambda )}.$$
(3.41b)
with $`w`$ defined as in (3.16). The second set of arguments that were added to $`A`$ and $`B`$ make explicit the position $`(\lambda _0,\mu _0)`$ of the source point which is causing the stresses at the field point $`(\lambda ,\mu )`$.

#### 3.2.6 The disc solution

The solution of equations (3.24) with right hand sides of the simplified form

$$\stackrel{~}{g}_1(\lambda ,\mu )=\delta (\lambda _0\lambda )\delta (\mu _0\mu ),\stackrel{~}{g}_2(\lambda ,\mu )=0,$$
(3.42)
is obtained from the solution (3.2.2) by interchanging $`\lambda \mu `$ and $`\lambda _0\mu _0`$. It is

$`S_{\lambda \lambda }`$ $`=`$ $`B(\mu ,\lambda ;\mu _0,\lambda _0)(\lambda _0\lambda )(\mu _0\mu )\delta (\mu _0\mu )(\lambda _0\lambda ),`$
$`S_{\mu \mu }`$ $`=`$ $`A(\mu ,\lambda ;\mu _0,\lambda _0)(\lambda _0\lambda )(\mu _0\mu ).`$
To find the solution to the full equations (3.24) at $`(\lambda ,\mu )`$, we multiply the singular solutions (3.2.2) and (3.2.6) by $`g_1(\lambda _0,\mu _0)`$ and $`g_2(\lambda _0,\mu _0)`$ respectively and integrate over $`D`$, the domain of dependence of $`(\lambda ,\mu )`$. This gives the first two lines of the two equations (3.44) below. The terms in the third lines are due to the boundary values of $`S_{\mu \mu }`$ at $`\mu =\alpha `$. They are found by multiplying the singular solution (3.2.2) evaluated for $`\mu _0=\alpha `$ by $`S_{\mu \mu }(\lambda _0,\alpha )`$ and integrating over $`\lambda _0`$ in $`D`$. It is easily verified that this procedure correctly represents the boundary values with singular solutions. The final result for the general solution of the Jeans equations (3.24) for Stäckel discs, after using the evaluations (3.41), is

$$\begin{array}{c}S_{\lambda \lambda }(\lambda ,\mu )=\underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0g_1(\lambda _0,\mu )\hfill \\ \hfill +\underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\mu }{\overset{\alpha }{}}d\mu _0\left[g_1(\lambda _0,\mu _0)\frac{2wE^{}(w)}{\pi (\mu _0\mu )}+g_2(\lambda _0,\mu _0)\frac{E(w)}{\pi (\lambda _0\mu )}\right]\\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0S_{\mu \mu }(\lambda _0,\alpha )[\frac{E(w)}{\pi (\lambda _0\mu )}\underset{\mu _0=\alpha }{]},\end{array}$$
(3.44a)
$$\begin{array}{c}S_{\mu \mu }(\lambda ,\mu )=\underset{\mu }{\overset{\alpha }{}}d\mu _0g_2(\lambda ,\mu _0)\hfill \\ \hfill +\underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\mu }{\overset{\alpha }{}}d\mu _0\left[g_1(\lambda _0,\mu _0)\frac{E(w)}{\pi (\lambda \mu _0)}g_2(\lambda _0,\mu _0)\frac{2wE^{}(w)}{\pi (\lambda _0\lambda )}\right]\\ \hfill +S_{\mu \mu }(\lambda ,\alpha )\underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0S_{\mu \mu }(\lambda _0,\alpha )[\frac{2wE^{}(w)}{\pi (\lambda _0\lambda )}\underset{\mu _0=\alpha }{]}.\end{array}$$
(3.44b)
The terms $`(\mu _0\mu )^1`$ and $`(\lambda _0\lambda )^1`$ do not cause singularities because they are canceled by components of $`w`$. In order to show that equations (3.44) are equivalent to the solution (3.21) given by Riemann’s method, integrate the terms in $`E^{}(w)`$ by parts, and use the definitions of $`S_{\tau \tau }`$, $`g_1`$ and $`g_2`$.

#### 3.2.7 Convergence of the disc solution

We now return to the convergence issues first discussed in §3.1.4, where we assumed that the density $`\rho `$ decays as $`N(\mu )\lambda ^{s/2}`$ at large distances and the Stäckel potential as $`𝒪(\lambda ^\delta )`$. For the physical reasons given there, the assigned boundary stress $`T_{\mu \mu }(\lambda ,\alpha )`$ cannot exceed $`𝒪(\lambda ^{\delta s/2})`$ at large $`\lambda `$, giving an $`S_{\mu \mu }(\lambda ,\alpha )`$ of $`𝒪(\lambda ^{\delta s/2+1/2})`$. It follows that the infinite integrals in $`S_{\mu \mu }(\lambda _0,\alpha )`$ in the solution (3.44) require only that $`s>2\delta +1`$ for their convergence. This is the less restrictive result to which we referred earlier.

The terms in the boundary stress are seen to contribute terms of the correct order $`𝒪(\lambda ^{\delta s/2+1/2})`$ to $`S_{\lambda \lambda }(\lambda ,\mu )`$ and $`S_{\mu \mu }(\lambda ,\mu )`$. The formulas for the density and potential show that $`g_1(\lambda ,\mu )=𝒪(\lambda ^{\delta s/21/2})`$ while $`g_2(\lambda ,\mu )`$ is larger and $`𝒪(\lambda ^{s/21/2})`$ as $`\lambda \mathrm{}`$. The $`\lambda _0`$ integrations with $`g_1`$ and $`g_2`$ in their integrands all converge provided $`s>2\delta +1`$. Hence, both $`S_{\lambda \lambda }(\lambda ,\mu )`$ and $`S_{\mu \mu }(\lambda ,\mu )`$ are $`𝒪(\lambda ^{\delta s/2+1/2})`$, so that the stress components $`T_{\tau \tau }(\lambda ,\mu )`$ ($`\tau =\lambda ,\mu `$) are $`𝒪(\lambda ^{\delta s/2})`$, which is consistent with the physical reasoning of §3.1.4.

Hence, all the conditions necessary for (3.44) to be a valid solution of the Jeans equations (3.24) for a Stäckel disc are satisfied provided that $`s>2\delta +1`$. We have seen in §3.1.4 that $`\delta `$ must lie in the range $`[\frac{1}{2},0)`$. When $`\delta 0`$ the models approach the isothermal disk, for which also $`s=1`$ when the density is consistent with the potential. Only then our requirement $`s>2\delta +1`$ is violated.

### 3.3 Alternative boundary conditions

We now derive the alternative form of the general disc solution when the boundary conditions are not specified on $`\mu =\alpha `$ but on $`\mu =\beta `$, or on $`\lambda =\alpha `$ rather than in the limit $`\lambda \mathrm{}`$. While the former switch is straightforward, the latter is non-trivial, and leads to non-physical solutions.

#### 3.3.1 Boundary condition for $`\mu `$

The analysis in §3.1 and §3.2 is that needed when the boundary conditions are imposed at large $`\lambda `$ and at $`\mu =\alpha `$. The Jeans equations (2.5.3) can be solved in a similar way when one or both of those conditions are imposed instead at the opposite boundaries $`\lambda =\alpha `$ and/or $`\mu =\beta `$. The solution by Riemann’s method is accomplished by applying Green’s theorem to a different domain, for example $`D^{}=\{(\lambda _0,\mu _0)`$: $`\lambda \lambda _0\mathrm{},\beta \mu _0\mu \}`$ when the boundary conditions are at $`\mu =\beta `$ and as $`\lambda \mathrm{}`$. The Riemann–Green functions have to satisfy the same PDE (3.10) and the same boundary conditions (3.12) and (3.1.1), and so again are given by equations (3.20a) and (3.20b). The variable $`w`$ is negative in $`D^{}`$ instead of positive as in $`D`$, but this is unimportant. The only significant difference in the solution of eq. (3.4) is that of a sign due to changes in the limits of the line integrals. The final result, in place of eq. (3.14), is

$$\begin{array}{c}T(\lambda ,\mu )=\underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\beta }{\overset{\mu }{}}d\mu _0𝒢(\lambda _0,\mu _0)U(\lambda _0,\mu _0)\hfill \\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0[(\frac{T}{\lambda _0}+\frac{c_2T}{\lambda _0\mu _0})𝒢\underset{\mu _0=\beta }{]}.\end{array}$$
(3.45)
To apply the method of singular solutions to solve for the stresses when the boundary stresses are specified at $`\mu =\beta `$ rather than at $`\mu =\alpha `$, we modify the singular solutions (3.2.2) by replacing the step-function $`(\mu _0\mu )`$ by $`(\mu \mu _0)`$ throughout. No other change is needed because both functions give $`\delta (\mu \mu _0)`$ on partial differentiation with respect to $`\mu `$. The two-dimensional problem for $`A`$ and $`B`$ remains the same, and so, as with Riemann’s method, its solution remains the same. Summing over sources in $`D^{}`$ now gives

$$\begin{array}{c}S_{\lambda \lambda }(\lambda ,\mu )=\underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0g_1(\lambda _0,\mu )\hfill \\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\beta }{\overset{\mu }{}}d\mu _0\left[g_1(\lambda _0,\mu _0)\frac{2wE^{}(w)}{\pi (\mu _0\mu )}+g_2(\lambda _0,\mu _0)\frac{E(w)}{\pi (\lambda _0\mu )}\right]\\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0S_{\mu \mu }(\lambda _0,\beta )[\frac{E(w)}{\pi (\lambda _0\mu )}\underset{\mu _0=\beta }{]},\end{array}$$
(3.46a)
$$\begin{array}{c}S_{\mu \mu }(\lambda ,\mu )=\underset{\beta }{\overset{\mu }{}}d\mu _0g_2(\lambda ,\mu _0)\hfill \\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\beta }{\overset{\mu }{}}d\mu _0\left[g_1(\lambda _0,\mu _0)\frac{E(w)}{\pi (\lambda \mu _0)}g_2(\lambda _0,\mu _0)\frac{2wE^{}(w)}{\pi (\lambda _0\lambda )}\right]\\ \hfill +S_{\mu \mu }(\lambda ,\beta )\underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0S_{\mu \mu }(\lambda _0,\beta )[\frac{2wE^{}(w)}{\pi (\lambda _0\lambda )}\underset{\mu _0=\beta }{]}.\end{array}$$
(3.46b)
as an alternative to equations (3.44).

#### 3.3.2 Boundary condition for $`\lambda `$

There is a much more significant difference when one assigns boundary values at $`\lambda =\alpha `$ rather than at $`\lambda \mathrm{}`$. It is still necessary that stresses decay to zero at large distances. The stresses induced by arbitrary boundary data at the finite boundary $`\lambda =\alpha `$ do decay to zero as a consequence of geometric divergence. The issue is that of the rate of this decay. We find that it is generally less than that required by our analysis in §3.1.4.

To isolate the effect of boundary data at $`\lambda =\alpha `$, we study solutions of the two-dimensional Jeans equations (2.5.3) when the inhomogeneous right hand side terms are set to zero and homogeneous boundary conditions of zero stress are applied at either $`\mu =\alpha `$ or $`\mu =\beta `$. These solutions can be derived either by Riemann’s method or by singular solutions. The solution of the homogeneous PDE $`T=0`$ is

$$T(\lambda ,\mu )=\underset{\mu }{\overset{\alpha }{}}\mathrm{d}\mu _0[(\frac{T}{\mu _0}\frac{c_1T}{\lambda _0\mu _0})𝒢(\lambda ,\mu ;\lambda _0,\mu _0)\underset{\lambda _0=\alpha }{]},$$
(3.47)
for the case of zero stress at $`\mu =\alpha `$, and

$$T(\lambda ,\mu )=\underset{\beta }{\overset{\mu }{}}\mathrm{d}\mu _0[(\frac{T}{\mu _0}\frac{c_1T}{\lambda _0\mu _0})𝒢(\lambda ,\mu ;\lambda _0,\mu _0)\underset{\lambda _0=\alpha }{]},$$
(3.48)
for the case of zero stress at $`\mu =\beta `$.

The behaviour of the stresses at large distances is governed by the behaviour of the Riemann–Green functions $`𝒢`$ for distant field points $`(\lambda ,\mu )`$ and source points at $`\lambda _0=\alpha `$. It follows from equations (3.20) that $`T_{\lambda \lambda }(\lambda ,\mu )=𝒪(\lambda ^{1/2})`$ and $`T_{\mu \mu }(\lambda ,\mu )=𝒪(\lambda ^{3/2})`$. As a restult, the radial stresses dominate at large distances and they decay as only the inverse first power of distance. Their rate of decay is less than $`𝒪(\lambda ^{\delta s/2})`$ – obtained in §3.1.4 from physical arguments – if the requirement $`s>2\delta +1`$ is satisfied. This inequality is the necessary condition which we derived in §3.2.6 for (3.44) to be a valid solution of the disc Jeans equations (3.24). It is violated in the isothermal limit.

There is a physical implication of radial stresses which decay as only the inverse first power of distance. It implies that net forces of finite magnitude are needed at an outer boundary to maintain the system, the finite magnitudes arising from the product of the decaying radial stresses and the increasing length of the boundary over which they act. That length grows as the first power of distance. Because this situation is perhaps more naturally understood in three dimensions, we return to it in our discussion of oblate models in §3.4.2. For now, lacking any physical reason for allowing a stellar system to have such an external constraint, we conclude that boundary conditions can be applied only at large $`\lambda `$ and not at $`\lambda =\alpha `$.

#### 3.3.3 Disc solution for a general finite region

We now apply the singular solution method to solve equations (3.24) in some rectangle $`\mu _{\mathrm{min}}\mu \mu _{\mathrm{max}}`$, $`\lambda _{\mathrm{min}}\lambda \lambda _{\mathrm{max}}`$, when the stress $`S_{\mu \mu }`$ is given a boundary in $`\mu `$, and $`S_{\lambda \lambda }`$ is given on a boundary in $`\lambda `$. This solution includes (3.44) and (3.46) as special cases. It will be needed for the large-radii scale-free case of §3.4.3.

As we saw in §3.3.1, singular solutions can easily be adapted to alternative choices for the domain of dependence of a field point $`(\lambda ,\mu )`$. Originally this was $`D`$, the first of the four quadrants into which $`(\lambda _0,\mu _0)`$-space is split by the lines $`\lambda _0=\lambda `$ and $`\mu _0=\mu `$ (Fig. 4). It has the singular solution (3.2.2). We then obtained the singular solution for the fourth quadrant $`D^{}`$ simply by replacing $`(\mu _0\mu )`$ by $`(\mu \mu _0)`$ in (3.2.2). We can similarly find the singular solution for the second quadrant $`\lambda _{\mathrm{min}}\lambda _0\lambda `$, $`\mu \mu _0\mu _{\mathrm{max}}`$ by replacing $`(\lambda _0\lambda )`$ by $`(\lambda \lambda _0)`$, and for the third quadrant $`\lambda _{\mathrm{min}}\lambda _0\lambda `$, $`\mu _{\mathrm{min}}\mu _0\mu `$ by replacing $`(\lambda _0\lambda )`$ by $`(\lambda \lambda _0)`$ and $`(\mu _0\mu )`$ by $`(\mu \mu _0)`$. We find the part of the solution of equations (3.24) due to the right hand side $`g`$ terms by multiplying the first and second terms of the singular solutions by $`g_1(\lambda _0,\mu _0)`$ and $`g_2(\lambda _0,\mu _0)`$, respectively, and integrating over the relevant domain. We use $`\lambda =\lambda _e`$ and $`\mu =\mu _e`$ to denote the boundaries at which stresses are specified. We find the part of the solution generated by the boundary values of $`S_{\mu \mu }`$ by multiplying the singular solution (3.2.2), modified for the domain and evaluated at $`\mu _0=\mu _e`$, by $`\pm S_{\mu \mu }(\lambda _0,\mu _e)`$ and integrating over $`\lambda _0`$ in the domain. The plus sign is needed when $`\mu _e=\mu _{\mathrm{min}}`$ and the minus when $`\mu _e=\mu _{\mathrm{max}}`$. Similarly, the part of the solution generated by the boundary values of $`S_{\lambda \lambda }`$ is obtained by multiplying the singular solution (3.2.6), modified for the domain and evaluated at $`\lambda _0=\lambda _e`$, by $`\pm S_{\lambda \lambda }(\lambda _e,\mu _0)`$ and integrating over $`\mu _0`$ in the domain. The sign is plus if $`\lambda _e=\lambda _{\mathrm{min}}`$ and minus if $`\lambda _e=\lambda _{\mathrm{max}}`$. The final solution is

$$\begin{array}{c}S_{\lambda \lambda }(\lambda ,\mu )=S_{\lambda \lambda }(\lambda _e,\mu )\underset{\lambda }{\overset{\lambda _e}{}}d\lambda _0g_1(\lambda _0,\mu )\hfill \\ \hfill +\underset{\lambda }{\overset{\lambda _e}{}}d\lambda _0\underset{\mu }{\overset{\mu _e}{}}d\mu _0\left[g_1(\lambda _0,\mu _0)B(\mu ,\lambda ;\mu _0,\lambda _0)+g_2(\lambda _0,\mu _0)A(\lambda ,\mu ;\lambda _0,\mu _0)\right]\\ \hfill \underset{\lambda }{\overset{\lambda _e}{}}d\lambda _0S_{\mu \mu }(\lambda _0,\mu _e)A(\lambda ,\mu ;\lambda _0,\mu _e)\underset{\mu }{\overset{\mu _e}{}}d\mu _0S_{\lambda \lambda }(\lambda _e,\mu _0)B(\mu ,\lambda ;\mu _0,\lambda _e),\end{array}$$
(3.49a)
$$\begin{array}{c}S_{\mu \mu }(\lambda ,\mu )=S_{\mu \mu }(\lambda ,\mu _e)\underset{\mu }{\overset{\mu _e}{}}d\mu _0g_2(\lambda ,\mu _0)\hfill \\ \hfill +\underset{\lambda }{\overset{\lambda _e}{}}d\lambda _0\underset{\mu }{\overset{\mu _e}{}}d\mu _0\left[g_1(\lambda _0,\mu _0)A(\mu ,\lambda ;\mu _0,\lambda _0)+g_2(\lambda _0,\mu _0)B(\lambda ,\mu ;\lambda _0,\mu _0)\right]\\ \hfill \underset{\lambda }{\overset{\lambda _e}{}}d\lambda _0S_{\mu \mu }(\lambda _0,\mu _e)B(\lambda ,\mu ;\lambda _0,\mu _e)\underset{\mu }{\overset{\mu _e}{}}d\mu _0S_{\lambda \lambda }(\lambda _e,\mu _0)A(\mu ,\lambda ;\mu _0,\lambda _e).\end{array}$$
(3.49b)
This solution is uniquely determined once $`g_1`$ and $`g_2`$ are given, and the boundary values $`S_{\mu \mu }(\lambda _0,\mu _e)`$ and $`S_{\lambda \lambda }(\lambda _e,\mu _0)`$ are prescribed. It shows that the hyperbolic equations (3.24) can equally well be integrated in either direction in the characteristic variables $`\lambda `$ and $`\mu `$. Solutions (3.44) and (3.46) are obtained by taking $`\lambda _e\mathrm{}`$, $`S_{\lambda \lambda }(\lambda _e,\mu _0)0`$, setting $`\mu _e=\alpha `$ and $`\mu _e=\beta `$ respectively, and evaluating $`A`$ and $`B`$ by equations (3.41).

### 3.4 Applying the disc solution to limiting cases

We showed in §2.6 that the Jeans equations for prolate and oblate potentials and for three-dimensional Stäckel models with a scale-free DF all reduce to a set of two equations equivalent to those for the Stäckel disc. Here we apply our solution for the Stäckel disc to these special three-dimensional cases, with particular attention to the behaviour at large radii and the boundary conditions. This provides further insight in some of the previously published solutions. We also consider the case of a Stäckel disc built with thin tube orbits.

#### 3.4.1 Prolate potentials

We can apply the disc solution (3.46) to solve the Jeans equations (2.6) by setting $`S_{\lambda \lambda }(\lambda ,\mu )=|\lambda \mu |^{\frac{1}{2}}𝒯_{\lambda \lambda }(\lambda ,\mu )`$ and $`S_{\mu \mu }(\lambda ,\mu )=|\mu \lambda |^{\frac{1}{2}}𝒯_{\mu \mu }(\lambda ,\mu )`$, and taking

$`g_1(\lambda ,\mu )`$ $`=`$ $`|\lambda \mu |^{\frac{1}{2}}(\lambda +\beta )^{\frac{1}{2}}(\mu +\beta )^{\frac{1}{2}}\left[\rho {\displaystyle \frac{V_S}{\lambda }}+{\displaystyle \frac{T_{\chi \chi }}{\lambda }}\right],`$
$`g_2(\lambda ,\mu )`$ $`=`$ $`|\mu \lambda |^{\frac{1}{2}}(\lambda +\beta )^{\frac{1}{2}}(\mu +\beta )^{\frac{1}{2}}\left[\rho {\displaystyle \frac{V_S}{\mu }}+{\displaystyle \frac{T_{\chi \chi }}{\mu }}\right].`$
The boundary terms in $`S_{\mu \mu }(\lambda ,\beta )`$ vanish because of the boundary condition (2.36). As before, we regard the azimuthal stress $`T_{\chi \chi }`$ as a variable that can be arbitrarily assigned, provided that it has the correct behaviour at large $`\lambda `$3.1.4). The choice of $`T_{\chi \chi }`$ is also restricted by the requirement that the resulting solutions for the stresses $`T_{\lambda \lambda }`$ and $`T_{\mu \mu }`$ must be non-negative (see §2.3).

The analysis needed to show that the solution obtained in this way is valid requires only minor modifications of that of §3.2.7. We suppose that the prescribed azimuthal stresses also decay as $`𝒪(\lambda ^{\delta s/2})`$ as $`\lambda \mathrm{}`$. As a result of the extra factor in the definitions (3.4.1), we now have $`g_1(\lambda ,\mu )=𝒪(\lambda ^{\delta s/2})`$ and $`g_2(\lambda ,\mu )=𝒪(\lambda ^{s/2})`$ as $`\lambda \mathrm{}`$. The $`\lambda _0`$ integrations converge provided $`s>2\delta +2`$, and $`S_{\lambda \lambda }`$ and $`S_{\mu \mu }`$ are $`𝒪(\lambda ^{\delta s/2+1})`$. Hence the stresses $`T_{\lambda \lambda }`$ and $`T_{\mu \mu }`$, which follow from $`T_{\tau \tau }=T_{\chi \chi }+S_{\tau \tau }/\sqrt{(\lambda \mu )(\lambda +\beta )(\mu +\beta )}`$, are once again $`𝒪(\lambda ^{\delta s/2})`$. The requirement $`s>2\delta +2`$ is no stronger than the requirement $`s>2\delta +1`$ of §3.2.7; it is simply the three-dimensional version of that requirement. It also does not break down until the isothermal limit. That limit is still $`\delta 0`$, but now $`s2`$.

#### 3.4.2 Oblate potentials

The oblate case with Jeans equations (2.6) differs significantly from the prolate case. Now $`S_{\lambda \lambda }(\lambda ,\nu )=|\lambda \nu |^{\frac{1}{2}}𝒯_{\lambda \lambda }(\lambda ,\nu )`$ vanishes at $`\lambda =\alpha `$ and $`S_{\nu \nu }(\lambda ,\nu )=|\nu \lambda |^{\frac{1}{2}}𝒯_{\nu \nu }(\lambda ,\nu )`$ vanishes at $`\nu =\alpha `$. If one again supposes that the azimuthal stresses $`T_{\varphi \varphi }`$ can be assigned initially, then one encounters the problem discussed in §3.3.2 of excessively large radial stresses at large distances. To relate that analysis to the present case, we use the solution (3.44) with $`\mu `$ replaced by $`\nu `$, and with zero boundary value $`S_{\nu \nu }(\lambda ,\alpha )`$, and for $`g_1`$ and $`g_2`$ the right hand side of (2.6) multiplied by $`|\lambda \nu |^{\frac{1}{2}}`$ and $`|\nu \lambda |^{\frac{1}{2}}`$, respectively.

The estimates we obtained for the prolate case are still valid, so the stresses $`T_{\lambda \lambda }`$ and $`T_{\nu \nu }`$ are $`𝒪(\lambda ^{\delta s/2})`$. Difficulties arise when this solution for $`S_{\lambda \lambda }`$ does not vanish at $`\lambda =\alpha `$, but instead has some nonzero value $`\kappa (\nu )`$ there. To obtain a physically acceptable solution, we must add to it a solution of the homogeneous equations (2.6) with boundary values $`𝒯_{\lambda \lambda }(\alpha ,\nu )=\kappa (\nu )/\sqrt{\alpha \nu }`$ and $`𝒯_{\nu \nu }(\lambda ,\alpha )=0`$. This is precisely the problem we discussed in §3.3.2 where we showed that the resulting solution gives $`𝒯_{\lambda \lambda }(\lambda ,\mu )=𝒪(\lambda ^{1/2})`$, and hence $`T_{\lambda \lambda }(\lambda ,\mu )=𝒪(\lambda ^1)`$. This is larger than $`𝒪(\lambda ^{\delta s/2})`$ when the three-dimensional requirement $`s>2\delta +2`$ is met. We therefore conclude that the approach in which one first selects the azimuthal stress $`T_{\varphi \varphi }`$ and then calculates the other two stresses will be unsuccessful unless the choice of $`T_{\varphi \varphi }`$ is fortunate, and leads to $`\kappa (\nu )0`$. Otherwise, it leads only to models which either violate the continuity condition $`T_{\lambda \lambda }T_{\varphi \varphi }=0`$ at $`\lambda =\alpha `$, or else have radial stresses which require external forces at large distances.

The physical implication of radial stresses which decay as only $`𝒪(\lambda ^1)`$, or the inverse second power of distance, is that net forces of finite magnitude are needed at an outer boundary to maintain the system. This finite magnitude arises from the product of the decaying radial stresses and the increasing surface area of the boundary over which they act, which grows as the second power of distance. This situation is analogous to that of an isothermal sphere, as illustrated in problem 4–9 of Binney & Tremaine (1987), for which the contribution from an outer surface integral must be taken into account in the balance between energies required by the virial theorem.

There are, of course, many physical models which satisfy the continuity condition and whose radial stresses decay in the physically correct manner at large distances, but some strategy other than that of assigning $`T_{\varphi \varphi }`$ initially is needed to find them. In fact, only Evans (1990) used the approach of assigning $`T_{\varphi \varphi }`$ initially. He computed a numerical solution for a mass model with $`s=3`$ and $`V_S𝒪(\lambda ^{1/2}\mathrm{ln}\lambda )`$ for large $`\lambda `$, so that the stresses there should be $`𝒪(\lambda ^2\mathrm{ln}\lambda )`$. He set $`T_{\varphi \varphi }=\frac{1}{3}\rho V_S`$, which is of this magnitude, and integrated from $`\lambda =\alpha `$ in the direction of increasing $`\lambda `$ for a finite range. Evans does not report on the large $`\lambda `$ behaviour, and it is possible that his choice of $`T_{\varphi \varphi }`$ gives $`\kappa (\nu )=0`$, but his Figure 2 especially shows velocity ellipsoids which become increasingly elongated in the radial direction, consistent with our prediction that $`T_{\lambda \lambda }`$ generally grows as $`𝒪(\lambda ^1)`$ when the boundary value of $`T_{\lambda \lambda }`$ is assigned at $`\lambda =\alpha `$.

A more common and effective approach to solve the Jeans equations for oblate models has been to specify the ratio $`T_{\lambda \lambda }/T_{\nu \nu }`$, and then to solve for one of those stresses and $`T_{\varphi \varphi }`$ (Bacon, Simien & Monnet 1983; Dejonghe & de Zeeuw 1988; Evans & Lynden–Bell 1991; Arnold 1995). This leads to a much simpler mathematical problem with just a single first-order PDE. The characteristics of that PDE have non-negative slopes $`d\lambda /d\nu `$, and therefore cut across the coordinate lines of constant $`\lambda `$ and $`\nu `$. The solution is obtained by integrating in along the characteristics from large $`\lambda `$. The continuity conditions (2.5.2) are taken care of automatically, the region $`\gamma \nu \alpha \mathrm{}`$ is covered, and it is easy to verify that the stresses so obtained are everywhere positive.

#### 3.4.3 Large radii limit with scale-free DF

We found in §2.5.4 that the first of the Jeans equations in conical coordinates (2.5.4) reduces to an algebraic relation for the radial stress $`T_{rr}`$. The problem that remains is that of solving the second and third Jeans equations for $`T_{\mu \mu }`$ and $`T_{\nu \nu }`$. Those equations are exactly the same as those of the disc case after we apply the coordinate permutation $`\lambda \mu \nu `$, and the physical domain is $`\gamma \nu \beta \mu \alpha `$ with finite ranges of both variables. Hence, the solution (3.49) can be applied with $`T_{\mu \mu }`$ assigned at either $`\mu _e=\alpha `$ or $`\mu _e=\beta `$, and $`T_{\nu \nu }`$ at either $`\nu _e=\beta `$ or $`\nu _e=\gamma `$. For $`g_1`$ and $`g_2`$ we take the same expressions as for the disc case, i.e., the right-hand side of (3.24), but with $`\lambda \mu \nu `$ and multiplied by $`r^\zeta `$. To obtain $`T_{\mu \mu }`$ and $`T_{\nu \nu }`$ from the $`S_{\lambda \lambda }`$ and $`S_{\mu \mu }`$ respectively, we use the transformation

$$S_{\tau \tau }=(\mu \nu )^{\frac{1}{2}}r^\zeta T_{\tau \tau },\tau =\mu ,\nu ,$$
(3.51)
with $`\zeta >0`$ the scaling factor. We can choose to specify the stress components on the two boundaries $`\mu =\beta `$ and $`\nu =\beta `$. For a given radius $`r`$ these boundaries cover the circular cross section with the $`(x,z)`$-plane (Fig. 3). We can consider the $`(x,z)`$-plane as the starting space for the solution. It turns out that the latter also applies to the triaxial solution (§4.6.3) and compares well with Schwarzschild (1993), who used the same plane to start his numerically calculated orbits from.

#### 3.4.4 Thin tube orbits

For infinitesimally thin tube orbits in Stäckel discs we have that $`S_{\lambda \lambda }0`$2.5.6), so that equations (3.24) reduce to

$$\frac{S_{\mu \mu }}{2(\lambda \mu )}=g_1(\lambda ,\mu ),\frac{S_{\mu \mu }}{\mu }=g_2(\lambda ,\mu ).$$
(3.52)
A solution is possible only if the right hand side terms satisfy the subsidiary equation

$$g_2(\lambda ,\mu )=2\frac{}{\mu }\left[(\lambda \mu )g_1(\lambda ,\mu )\right].$$
(3.53)
We find below that this equation places restrictions on the form of the (surface) density $`\rho `$, and we use this relation between $`g_1`$ and $`g_2`$ to show that the disc solution (3.44) yields the right results for the stress components.

If we write the disc potential (2.24) as a divided difference, $`V_S=f[\lambda ,\mu ]`$, we have that

$$g_1=(\lambda \mu )^{\frac{1}{2}}\rho f[\lambda ,\lambda ,\mu ],g_2=(\lambda \mu )^{\frac{1}{2}}\rho f[\lambda ,\mu ,\mu ].$$
(3.54)
Upon substitution of these expressions in (3.53) we obtain a PDE in $`\mu `$, of which the solution implies the following form for the density

$$\rho (\lambda ,\mu )=\frac{\stackrel{~}{f}(\lambda )}{(\lambda \mu )\sqrt{f[\lambda ,\lambda ,\mu ]}},$$
(3.55)
where $`\stackrel{~}{f}(\lambda )`$ is an arbitrary function independent of $`\mu `$. From (3.52) and the definition (3.2.1) it then follows that $`T_{\mu \mu }(\lambda ,\mu ,\nu )=2\stackrel{~}{f}(\lambda )\sqrt{f[\lambda ,\lambda ,\mu ]}`$. The tube density that de Zeeuw, Hunter & Schwarzschild (1987) derive from the DF for thin tube orbits in the perfect elliptic disk (their eq. \[4.25\]) is indeed of the form (3.55).

To show that the general disc solution (3.44) gives $`S_{\lambda \lambda }(\lambda ,\mu )=0`$, we substitute eq. (3.53) for $`g_2(\lambda ,\mu )`$ in (3.44a). After partial integration and using

$$2(\lambda _0\mu _0)\frac{}{\mu _0}\frac{E(w)}{\pi (\lambda _0\mu )}=\frac{2wE^{}(w)}{\pi (\mu _0\mu )},$$
(3.56)
we find that the area integral reduces to

$$\underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0\{g_1(\lambda _0,\mu )2(\lambda _0+\alpha )g_1(\lambda _0,\alpha )\left[\frac{E(w)}{\pi (\lambda _0\mu )}\underset{\mu _0=\alpha }{]}\right\}.$$
(3.57)
The first part cancels the first line of (3.44a) and since from (3.52) we have that $`2(\lambda _0+\alpha )g_1(\lambda _0,\alpha )=S_{\mu \mu }(\lambda _0,\alpha )`$, the second part cancels the third line. Hence, we have $`S_{\lambda \lambda }(\lambda ,\mu )=0`$ as required. To see that the general disc solution also yields $`S_{\mu \mu }(\lambda ,\mu )`$ correctly, we apply similar steps to (3.44b), where we use the relation

$$2(\lambda _0\mu _0)\frac{}{\mu _0}\frac{2wE^{}(w)}{\pi (\lambda _0\lambda )}=\frac{E(w)}{\pi (\lambda \mu _0)}.$$
(3.58)
We are finally left with

$$S_{\mu \mu }(\lambda ,\mu )=S_{\mu \mu }(\lambda ,\alpha )\underset{\mu }{\overset{\alpha }{}}d\mu _0g_2(\lambda ,\mu _0),$$
(3.59)
which is just the second equation of (3.52) integrated with respect to $`\mu `$.

## 4 The general case

We now solve the system of three Jeans equations (2.16) for triaxial Stäckel models by applying the singular solution superposition method, introduced in §3.2 for the two-dimensional case. Although the calculations are more complex for a triaxial model, the step-wise solution method is similar to that in two dimensions. Specifically, we first simplify the Jeans equations and show that they reduce to a three-dimensional homogeneous boundary problem. We then find a two-parameter particular solution and apply contour integration to both complex parameters to obtain the general homogeneous solution. The latter yields the three singular solutions of the simplified Jeans equations, from which, by superposition, we construct the general solution.

### 4.1 Simplified Jeans equations

We start by introducing the functions

$$S_{\tau \tau }(\lambda ,\mu ,\nu )=\sqrt{(\lambda \mu )(\lambda \nu )(\mu \nu )}T_{\tau \tau }(\lambda ,\mu ,\nu ),$$
(4.1)
with $`\tau =\lambda ,\mu ,\nu `$, to write the Jeans equations for triaxial Stäckel models (2.16) in the more convenient form

$$\frac{S_{\lambda \lambda }}{\lambda }\frac{S_{\mu \mu }}{2(\lambda \mu )}\frac{S_{\nu \nu }}{2(\lambda \nu )}=g_1(\lambda ,\mu ,\nu ),$$
(4.2a)
$$\frac{S_{\mu \mu }}{\mu }\frac{S_{\nu \nu }}{2(\mu \nu )}\frac{S_{\lambda \lambda }}{2(\mu \lambda )}=g_2(\lambda ,\mu ,\nu ),$$
(4.2b)
$$\frac{S_{\nu \nu }}{\nu }\frac{S_{\lambda \lambda }}{2(\nu \lambda )}\frac{S_{\mu \mu }}{2(\nu \mu )}=g_3(\lambda ,\mu ,\nu ),$$
(4.2c)
where the function $`g_1`$ is defined as

$$g_1(\lambda ,\mu ,\nu )=\sqrt{(\lambda \mu )(\lambda \nu )(\mu \nu )}\rho \frac{V_S}{\lambda },$$
(4.3)
and $`g_2`$ and $`g_3`$ follow by cyclic permutation $`\lambda \mu \nu \lambda `$. We keep the three terms $`\lambda \mu `$, $`\lambda \nu `$ and $`\mu \nu `$ under one square root. With each cyclic permutation two of the three terms change sign, so that the combination of the three terms is always positive real. Therefore. the square root of the combination is always single-valued, whereas in the case of three separate square roots we would have a multi-valued function.

We simplify equations (4.2) by substituting for $`g_1`$, $`g_2`$ and $`g_3`$, respectively

$`\stackrel{~}{g}_1(\lambda ,\mu ,\nu )`$ $`=`$ $`0,`$
$`\stackrel{~}{g}_2(\lambda ,\mu ,\nu )`$ $`=`$ $`\delta (\lambda _0\lambda )\delta (\mu _0\mu )\delta (\nu _0\nu ),`$ (4.4)
$`\stackrel{~}{g}_3(\lambda ,\mu ,\nu )`$ $`=`$ $`0,`$
with

$$\gamma \nu \nu _0\beta \mu \mu _0\alpha \lambda \lambda _0.$$
(4.5)
We obtain two similar systems of simplified equations by cyclic permutation of the left-hand side of (4.2). Once we have obtained the singular solutions of the simplified system with the right-hand side given by (4.1), those for the other two systems follow via cyclic permutation.

### 4.2 Homogeneous boundary problem

The choice (4.1) implies that the functions $`S_{\tau \tau }(\lambda ,\mu ,\nu )`$ (4.1) must have the following forms

$`S_{\lambda \lambda }`$ $`=`$ $`A(\lambda ,\mu ,\nu )(\lambda _0\lambda )(\mu _0\mu )(\nu _0\nu )`$
$`+F(\lambda ,\mu )\delta (\nu _0\nu )(\lambda _0\lambda )(\mu _0\mu ),`$
$`S_{\mu \mu }`$ $`=`$ $`B(\lambda ,\mu ,\nu )(\lambda _0\lambda )(\mu _0\mu )(\nu _0\nu )`$
$`+G(\lambda ,\mu )\delta (\nu _0\nu )(\lambda _0\lambda )(\mu _0\mu )`$
$`+H(\mu ,\nu )\delta (\lambda _0\lambda )(\mu _0\mu )(\nu _0\nu )`$
$`\delta (\lambda _0\lambda )\delta (\nu _0\nu )(\mu _0\mu ),`$
$`S_{\nu \nu }`$ $`=`$ $`C(\lambda ,\mu ,\nu )(\lambda _0\lambda )(\mu _0\mu )(\nu _0\nu )`$
$`+I(\mu ,\nu )\delta (\lambda _0\lambda )(\mu _0\mu )(\nu _0\nu ),`$
with $`A`$, $`B`$, $`C`$ and $`F`$, $`G`$, $`H`$, $`I`$ yet unknown functions of three and two coordinates, respectively, and $``$ the step-function (3.26). After substituting these forms into the simplified Jeans equations and matching terms we obtain 14 equations. Eight of them comprise the following two homogeneous systems with two boundary conditions each

$$\{\begin{array}{ccccccc}\frac{F}{\lambda }\frac{G}{2(\lambda \mu )}\hfill & =& 0,\hfill & & F(\lambda _0,\mu )\hfill & =& \frac{1}{2(\lambda _0\mu )},\hfill \\ \frac{G}{\mu }\frac{F}{2(\mu \lambda )}\hfill & =& 0,\hfill & & G(\lambda ,\mu _0)\hfill & =& 0,\hfill \end{array}$$
(4.7)
and

$$\{\begin{array}{ccccccc}\frac{H}{\mu }\frac{I}{2(\mu \nu )}\hfill & =& 0,\hfill & & H(\mu _0,\nu )\hfill & =& 0,\hfill \\ \frac{I}{\nu }\frac{H}{2(\nu \mu )}\hfill & =& 0,\hfill & & I(\mu ,\nu _0)\hfill & =& \frac{1}{2(\nu _0\mu )}.\hfill \end{array}$$
(4.8)
We have shown in §3 how to solve these two-dimensional homogeneous boundary problems in terms of the complete elliptic integral of the second kind $`E`$ and its derivative $`E^{}`$. The solutions are

$`F(\lambda ,\mu )`$ $`={\displaystyle \frac{E(w)}{\pi (\lambda _0\mu )}},`$ $`G(\lambda ,\mu )`$ $`={\displaystyle \frac{2wE^{}(w)}{\pi (\lambda _0\lambda )}},`$
$`H(\mu ,\nu )`$ $`={\displaystyle \frac{2uE^{}(u)}{\pi (\nu _0\nu )}},`$ $`I(\mu ,\nu )`$ $`={\displaystyle \frac{E(u)}{\pi (\mu \nu _0)}},`$
where $`u`$ and similarly $`v`$, which we will encounter later on, follow from $`w`$ (3.16) by cyclic permutation $`\lambda \mu \nu \lambda `$ and $`\lambda _0\mu _0\nu _0\lambda _0`$, so that

$$u=\frac{(\mu _0\mu )(\nu _0\nu )}{(\mu _0\nu _0)(\mu \nu )},v=\frac{(\nu _0\nu )(\lambda _0\lambda )}{(\lambda _0\nu _0)(\lambda \nu )}.$$
(4.10)
The remaining six equations form a three-dimensional homogeneous boundary problem, consisting of three homogeneous Jeans equations

$`{\displaystyle \frac{A}{\lambda }}{\displaystyle \frac{B}{2(\lambda \mu )}}{\displaystyle \frac{C}{2(\lambda \nu )}}`$ $`=`$ $`0,`$
$`{\displaystyle \frac{B}{\mu }}{\displaystyle \frac{C}{2(\mu \nu )}}{\displaystyle \frac{A}{2(\mu \lambda )}}`$ $`=`$ $`0,`$ (4.11)
$`{\displaystyle \frac{C}{\nu }}{\displaystyle \frac{A}{2(\nu \lambda )}}{\displaystyle \frac{B}{2(\nu \mu )}}`$ $`=`$ $`0.`$
and three boundary conditions, specifically the values of $`A(\lambda _0,\mu ,\nu )`$, $`B(\lambda ,\mu _0,\nu )`$, and $`C(\lambda ,\mu ,\nu _0)`$. As in §3.2.2, it is useful to supplement these boundary conditions with the values of $`A`$, $`B`$, and $`C`$ at the other boundary surfaces. These are obtained by integrating the pairs of equations (4.2) which apply at those surfaces, and using the boundary conditions. This results in the following nine boundary values

$`A(\lambda _0,\mu ,\nu )`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}[{\displaystyle \frac{E(u)}{(\lambda _0\nu )(\mu \nu _0)}}+{\displaystyle \frac{2uE^{}(u)}{(\lambda _0\mu )(\nu _0\nu )}}],`$
$`A(\lambda ,\mu _0,\nu )`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}[{\displaystyle \frac{E(v)}{(\lambda _0\nu )(\mu _0\nu _0)}}+{\displaystyle \frac{2vE^{}(v)}{(\lambda _0\mu _0)(\nu _0\nu )}}],`$
$`A(\lambda ,\mu ,\nu _0)`$ $`=`$ $`{\displaystyle \frac{E(w)}{4\pi (\lambda _0\mu )}}[{\displaystyle \frac{\lambda \mu }{(\lambda \nu _0)(\mu \nu _0)}}+{\displaystyle \frac{\lambda _0\mu _0}{(\lambda _0\nu _0)(\mu _0\nu _0)}}],`$
$`B(\lambda _0,\mu ,\nu )`$ $`=`$ $`{\displaystyle \frac{uE^{}(u)}{2\pi (\nu _0\nu )}}[{\displaystyle \frac{\mu _0\mu }{(\lambda _0\mu _0)(\lambda _0\mu )}}{\displaystyle \frac{\nu _0\nu }{(\lambda _0\nu _0)(\lambda _0\nu )}}],`$
$`B(\lambda ,\mu _0,\nu )`$ $`=`$ $`\mathrm{\hspace{0.17em}0},`$ (4.12)
$`B(\lambda ,\mu ,\nu _0)`$ $`=`$ $`{\displaystyle \frac{wE^{}(w)}{2\pi (\lambda _0\lambda )}}[{\displaystyle \frac{\mu _0\mu }{(\mu _0\nu _0)(\mu \nu _0)}}{\displaystyle \frac{\lambda _0\lambda }{(\lambda _0\nu _0)(\lambda \nu _0)}}],`$
$`C(\lambda _0,\mu ,\nu )`$ $`=`$ $`{\displaystyle \frac{E(u)}{4\pi (\mu \nu _0)}}[{\displaystyle \frac{\mu \nu }{(\lambda _0\mu )(\lambda _0\nu )}}+{\displaystyle \frac{\mu _0\nu _0}{(\lambda _0\mu _0)(\lambda _0\nu _0)}}],`$
$`C(\lambda ,\mu _0,\nu )`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}[{\displaystyle \frac{E(v)}{(\lambda _0\mu _0)(\lambda \nu _0)}}{\displaystyle \frac{2vE^{}(v)}{(\mu _0\nu _0)(\lambda _0\lambda )}}],`$
$`C(\lambda ,\mu ,\nu _0)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}[{\displaystyle \frac{E(w)}{(\lambda _0\mu )(\lambda \nu _0)}}{\displaystyle \frac{2wE^{}(w)}{(\mu \nu _0)(\lambda _0\lambda )}}].`$
If we can solve the three homogeneous equations (4.2) and satisfy the nine boundary expressions (4.2) simultaneously, then we obtain the singular solutions (4.2). By superposition, we can then construct the solution of the Jeans equations for triaxial Stäckel models.

### 4.3 Particular solution

By analogy with the two-dimensional case, we look for particular solutions of the homogeneous equations (4.2) and by superposition of these particular solutions we try to satisfy the boundary expressions (4.2) simultaneously, in order to obtain the homogeneous solution for $`A`$, $`B`$ and $`C`$.

#### 4.3.1 One-parameter particular solution

By substitution one can verify that

$$A^P(\lambda ,\mu ,\nu )=\frac{\sqrt{(\lambda \mu )(\lambda \nu )(\mu \nu )}}{(\lambda \mu )(\lambda \nu )}\frac{(z\lambda )}{(z\mu )(z\nu )},$$
(4.13)
with $`B^P`$ and $`C^P`$ following from $`A^P`$ by cyclic permutation, solves the homogeneous equations (4.2). To satisfy the nine boundary expressions (4.2), we could integrate this particular solution over its free parameter $`z`$, in the complex plane. From §3.2.4, it follows that, at the boundaries, this results in simple polynomials in $`(\lambda ,\mu ,\nu )`$ and $`(\lambda _0,\mu _0,\nu _0)`$. This means that the nine boundary expressions (4.2) cannot be satisfied, since in addition to these simple polynomials they also contain $`E`$ and $`E^{}`$. The latter are functions of one variable, so that at least one extra freedom is necessary. Hence, we look for a particular solution with two free parameters.

#### 4.3.2 Two-parameter particular solution

A particular solution with two free parameters $`z_1`$ and $`z_2`$ can be found by splitting the $`z`$-dependent terms of the one-parameter solution (4.13) into two similar parts and then relabelling them. The result is the following two-parameter particular solution

$`A^P`$ $`=`$ $`{\displaystyle \frac{\sqrt{(\lambda \mu )(\lambda \nu )(\mu \nu )}}{(\lambda \mu )(\lambda \nu )}}{\displaystyle \underset{i=1}{\overset{2}{}}}{\displaystyle \frac{(z_i\lambda )^{\frac{1}{2}}}{(z_i\mu )^{\frac{1}{2}}(z_i\nu )^{\frac{1}{2}}}},`$
$`B^P`$ $`=`$ $`{\displaystyle \frac{\sqrt{(\lambda \mu )(\lambda \nu )(\mu \nu )}}{(\mu \nu )(\mu \lambda )}}{\displaystyle \underset{i=1}{\overset{2}{}}}{\displaystyle \frac{(z_i\mu )^{\frac{1}{2}}}{(z_i\nu )^{\frac{1}{2}}(z_i\lambda )^{\frac{1}{2}}}},`$ (4.14)
$`C^P`$ $`=`$ $`{\displaystyle \frac{\sqrt{(\lambda \mu )(\lambda \nu )(\mu \nu )}}{(\nu \lambda )(\nu \mu )}}{\displaystyle \underset{i=1}{\overset{2}{}}}{\displaystyle \frac{(z_i\nu )^{\frac{1}{2}}}{(z_i\lambda )^{\frac{1}{2}}(z_i\mu )^{\frac{1}{2}}}}.`$
These functions are cyclic in $`(\lambda ,\mu ,\nu )`$, as is required from the symmetry of the homogeneous equations (4.2). The presence of the square roots, such as occurred earlier in the solution (3.2.3) for the disc case, allows us to fit boundary values that contain elliptic integrals.

To show that this particular solution solves the homogeneous Jeans equations, we calculate the derivative of $`A^P(\lambda ,\mu ,\nu )`$ with respect to $`\lambda `$:

$$\frac{A^P}{\lambda }=\frac{A^P}{2}\left(\frac{1}{\lambda z_1}+\frac{1}{\lambda z_2}\frac{1}{\lambda \mu }\frac{1}{\lambda \nu }\right).$$
(4.15)
This can be written as

$`{\displaystyle \frac{A^P}{\lambda }}`$ $`=`$ $`{\displaystyle \frac{1}{2(\lambda \mu )}}\left[{\displaystyle \frac{(z_1\mu )(z_2\mu )(\lambda \nu )}{(z_1\lambda )(z_2\lambda )(\mu \nu )}}A^P\right]`$
$`+{\displaystyle \frac{1}{2(\lambda \nu )}}\left[{\displaystyle \frac{(z_1\nu )(z_2\nu )(\lambda \mu )}{(z_1\lambda )(z_2\lambda )(\mu \nu )}}A^P\right].`$
From the two-parameter particular solution we have

$`{\displaystyle \frac{B^P}{A^P}}`$ $`=`$ $`{\displaystyle \frac{(z_1\mu )(z_2\mu )(\lambda \nu )}{(z_1\lambda )(z_2\lambda )(\mu \nu )}},`$
$`{\displaystyle \frac{C^P}{A^P}}`$ $`=`$ $`{\displaystyle \frac{(z_1\nu )(z_2\nu )(\lambda \mu )}{(z_1\lambda )(z_2\lambda )(\mu \nu )}},`$
so that, after substitution of these ratios, the first homogeneous equation of (4.2), is indeed satisfied. The remaining two homogeneous equations can be checked in the same way.

### 4.4 The homogeneous solution

In order to satisfy the four boundary expressions of the two-dimensional case, we multiplied the one-parameter particular solution by terms depending on $`\lambda _0`$, $`\mu _0`$ and the free complex parameter $`z`$, followed by contour integration over the latter. Similarly, in the triaxial case we multiply the two-parameter particular solution (3.35) by terms depending on $`\lambda _0`$, $`\mu _0`$, $`\nu _0`$ and the two free parameters $`z_1`$ and $`z_2`$, in such a way that by contour integration over the latter two complex parameters the nine boundary expressions (4.2) can be satisfied. Since these terms and the integration are independent of $`\lambda `$, $`\mu `$ and $`\nu `$, it follows from the superposition principle that the homogeneous equations (4.2) remain satisfied.

The contour integrations over $`z_1`$ and $`z_2`$ are mutually independent, since we can separate the two-parameter particular solution (4.3.2) with respect to these two parameters. This allows us to choose a pair of contours, one contour in the $`z_1`$-plane and the other contour in the $`z_2`$-plane, and integrate over them separately. We consider the same simple contours as in the disk case (Fig. 5) around the pairs of branch points $`(\lambda ,\lambda _0)`$ and $`(\mu ,\mu _0)`$, and a similar contour around $`(\nu ,\nu _0)`$. We denote these contours by $`C_i^\lambda `$, $`C_i^\mu `$ and $`C_i^\nu `$ respectively, with $`i=1,2`$ indicating in which of the two complex planes we apply the contour integration.

#### 4.4.1 Boundary expressions for $`B`$

It follows from (4.2) that $`B=0`$ at the boundary $`\mu =\mu _0`$. From Cauchy’s theorem, $`B`$ would indeed vanish if, in this case, in either the $`z_1`$-plane or $`z_2`$-plane the integrand for $`B`$ is analytic within the chosen integration contour. The boundary expression for $`B`$ at $`\nu =\nu _0`$ follows from the one at $`\lambda =\lambda _0`$ by taking $`\nu \lambda `$ and $`\nu _0\lambda _0`$. In addition to this symmetry, also the form of both boundary expressions puts constraints on the solution for $`B`$. The boundary expressions can be separated in two parts, one involving the complete elliptic integral $`E^{}`$ and the other consisting of a two-component polynomial in $`\tau `$ and $`\tau _0`$ ($`\tau =\lambda ,\mu ,\nu `$). Each of the two parts follows from a contour integration in one of the two complex planes. For either of the complex parameters, $`z_1`$ or $`z_2`$, the integrands will consist of a combination of the six terms $`z_i\tau `$ and $`z_i\tau _0`$ with powers that are half-odd integers, i.e., the integrals are of hyperelliptic form. If two of the six terms cancel on one of the boundaries, we will be left with an elliptic integral. We expect the polynomial to result from applying the Residue theorem to a double pole, as this would involve a first derivative and hence give two components. This leads to the following Ansatz

$$\begin{array}{c}B(\lambda ,\mu ,\nu )\frac{\sqrt{(\lambda \mu )(\lambda \nu )(\mu \nu )}}{(\mu \nu )(\mu \lambda )}\times \hfill \\ \hfill \underset{C_1}{}\frac{(z_1\mu )^{\frac{1}{2}}(z_1\lambda _0)^{\frac{1}{2}}\mathrm{d}z_1}{(z_1\nu )^{\frac{1}{2}}(z_1\lambda )^{\frac{1}{2}}(z_1\mu _0)^{\frac{1}{2}}(z_1\nu _0)^{\frac{3}{2}}}\times \\ \hfill \underset{C_2}{}\frac{(z_2\mu )^{\frac{1}{2}}(z_2\nu _0)^{\frac{1}{2}}\mathrm{d}z_2}{(z_2\nu )^{\frac{1}{2}}(z_2\lambda )^{\frac{1}{2}}(z_2\mu _0)^{\frac{1}{2}}(z_2\lambda _0)^{\frac{3}{2}}}.\end{array}$$
(4.18)
Upon substitution of $`\mu =\mu _0`$, the terms involving $`\mu _0`$ cancel in both integrals, so that the integrands are analytic in both contours $`C_1^\mu `$ and $`C_2^\mu `$. Hence, by choosing either of these contours as integration contour, the boundary expression $`B(\lambda ,\mu _0,\nu )=0`$ is satisfied.

When $`\lambda =\lambda _0`$, the terms with $`\lambda _0`$ in the first integral in (4.18) cancel, while in the second integral we have $`(z_2\lambda _0)^2`$. The first integral is analytic within $`C_1^\lambda `$, so that there is no contribution from this contour. However, the integral over $`C_1^\mu `$ is elliptic and can be evaluated in terms of $`E^{}`$ (cf. §3.2.5). We apply the Residue theorem to the second integral, for which there is a double pole inside the contour $`C_2^\lambda `$. Considering $`C_1^\mu `$ and $`C_2^\lambda `$ as a pair of contours, the expression for $`B`$ at $`\lambda =\lambda _0`$ becomes

$$\begin{array}{c}B(\lambda ,\mu ,\nu )16\pi ^2\frac{\sqrt{(\lambda _0\mu _0)(\lambda _0\nu _0)(\mu _0\nu _0)}}{(\mu _0\nu _0)(\mu _0\lambda _0)}\times \hfill \\ \hfill \frac{uE^{}(u)}{2\pi (\nu _0\nu )}[\frac{\mu _0\mu }{(\lambda _0\mu _0)(\lambda _0\mu )}\frac{\nu _0\nu }{(\lambda _0\nu _0)(\lambda _0\nu )}],\end{array}$$
(4.19)
which is the required boundary expression up to a scaling factor. As before, we keep the terms $`\lambda _0\mu _0`$, $`\lambda _0\nu _0`$ and $`\mu _0\nu _0`$ under one square root, so that it is single-valued with respect to cyclic permutation in these coordinates.

The boundary expression for $`B`$ at $`\nu =\nu _0`$ is symmetric with the one at $`\lambda =\lambda _0`$, so that a similar approach can be used. In this case, for the second integral, there is no contribution from $`C_2^\nu `$, whereas it can be expressed in terms of $`E^{}`$ if $`C_2=C_2^\mu `$. The first integrand has a double pole in $`C_1^\nu `$. The total contribution from the pair ($`C_1^\nu `$,$`C_2^\mu `$) gives the correct boundary expression, up to a scaling factor that is the same as in (4.19).

Taking into account the latter scaling factor, this shows that the Ansatz (4.18) for $`B`$ produces the correct boundary expressions and hence we postulate it as the homogeneous solution for $`B`$. The expressions for $`A`$ and $`C`$ then follow from the ratios (4.3.2). Absorbing the minus sign in (4.19) into the pair of contours, i.e., either of the two contours we integrate in clockwise direction, we postulate the following homogeneous solution

$`A(\lambda ,\mu ,\nu )`$ $`=`$ $`{\displaystyle \frac{(\mu _0\nu _0)(\mu _0\lambda _0)}{16\pi ^2(\lambda \mu )(\lambda \nu )}}\sqrt{{\displaystyle \frac{(\lambda \mu )(\lambda \nu )(\mu \nu )}{(\lambda _0\mu _0)(\lambda _0\nu _0)(\mu _0\nu _0)}}}\times `$
$`{\displaystyle \underset{C_1}{}}`$ $`{\displaystyle \frac{(z_1\lambda )^{\frac{1}{2}}(z_1\lambda _0)^{\frac{1}{2}}\mathrm{d}z_1}{(z_1\mu )^{\frac{1}{2}}(z_1\nu )^{\frac{1}{2}}(z_1\mu _0)^{\frac{1}{2}}(z_1\nu _0)^{\frac{3}{2}}}}\times `$
$`{\displaystyle \underset{C_2}{}}`$ $`{\displaystyle \frac{(z_2\lambda )^{\frac{1}{2}}(z_2\nu _0)^{\frac{1}{2}}\mathrm{d}z_2}{(z_2\mu )^{\frac{1}{2}}(z_2\nu )^{\frac{1}{2}}(z_2\mu _0)^{\frac{1}{2}}(z_2\lambda _0)^{\frac{3}{2}}}},`$ (4.20)
$`B(\lambda ,\mu ,\nu )`$ $`=`$ $`{\displaystyle \frac{(\mu _0\nu _0)(\mu _0\lambda _0)}{16\pi ^2(\mu \nu )(\mu \lambda )}}\sqrt{{\displaystyle \frac{(\lambda \mu )(\lambda \nu )(\mu \nu )}{(\lambda _0\mu _0)(\lambda _0\nu _0)(\mu _0\nu _0)}}}\times `$
$`{\displaystyle \underset{C_1}{}}`$ $`{\displaystyle \frac{(z_1\mu )^{\frac{1}{2}}(z_1\lambda _0)^{\frac{1}{2}}\mathrm{d}z_1}{(z_1\nu )^{\frac{1}{2}}(z_1\lambda )^{\frac{1}{2}}(z_1\mu _0)^{\frac{1}{2}}(z_1\nu _0)^{\frac{3}{2}}}}\times `$
$`{\displaystyle \underset{C_2}{}}`$ $`{\displaystyle \frac{(z_2\mu )^{\frac{1}{2}}(z_2\nu _0)^{\frac{1}{2}}\mathrm{d}z_2}{(z_2\nu )^{\frac{1}{2}}(z_2\lambda )^{\frac{1}{2}}(z_2\mu _0)^{\frac{1}{2}}(z_2\lambda _0)^{\frac{3}{2}}}},`$ (4.21)
$`C(\lambda ,\mu ,\nu )`$ $`=`$ $`{\displaystyle \frac{(\mu _0\nu _0)(\mu _0\lambda _0)}{16\pi ^2(\nu \lambda )(\nu \mu )}}\sqrt{{\displaystyle \frac{(\lambda \mu )(\lambda \nu )(\mu \nu )}{(\lambda _0\mu _0)(\lambda _0\nu _0)(\mu _0\nu _0)}}}\times `$
$`{\displaystyle \underset{C_1}{}}`$ $`{\displaystyle \frac{(z_1\nu )^{\frac{1}{2}}(z_1\lambda _0)^{\frac{1}{2}}\mathrm{d}z_1}{(z_1\lambda )^{\frac{1}{2}}(z_1\mu )^{\frac{1}{2}}(z_1\mu _0)^{\frac{1}{2}}(z_1\nu _0)^{\frac{3}{2}}}}\times `$
$`{\displaystyle \underset{C_2}{}}`$ $`{\displaystyle \frac{(z_2\nu )^{\frac{1}{2}}(z_2\nu _0)^{\frac{1}{2}}\mathrm{d}z_2}{(z_2\lambda )^{\frac{1}{2}}(z_2\mu )^{\frac{1}{2}}(z_2\mu _0)^{\frac{1}{2}}(z_2\lambda _0)^{\frac{3}{2}}}}.`$ (4.22)
The integrands consist of multi-valued functions that all come in pairs of the form $`(z\tau )^{\frac{1}{2}m}(z\tau _0)^{\frac{1}{2}n}`$, for integers $`m`$ and $`n`$, with $`\tau `$ equal to $`\lambda `$, $`\mu `$ or $`\nu `$. Hence, completely analogous to our procedure in §3.2.4, we can make the integrands single-valued by specifying, in the complex $`z_1`$-plane and $`z_2`$-plane, three cuts running between the three pairs $`(\lambda ,\lambda _0)`$, $`(\mu ,\mu _0)`$, $`(\nu ,\nu _0)`$ of branch points, that are enclosed by the simple contours. The integrands are now analytic in the cut plane away from its cuts and behave again as $`z_i^2`$ at large distances, so that the integral over a circular contour with radius going to infinity, will be zero. Hence, connecting the simple contours $`C_i^\lambda `$, $`C_i^\mu `$ and $`C_i^\nu `$ with this circular contour, shows that their cumulative contribution cancels

$$C_i^\nu +C_i^\mu +C_i^\lambda =0,i=1,2.$$
(4.23)
This relation will allow us to make a combination of contours, so that the nine boundary expressions (4.2) can be satisfied simultaneously4.4.3). Before doing so, we first establish whether, with the homogeneous solution for $`A`$ and $`C`$ given by (4.4.1) and (4.4.1), respectively, we indeed satisfy their corresponding boundary expressions separately.

#### 4.4.2 Boundary expressions for $`A`$ and $`C`$

The boundary expressions and the homogeneous solution of $`C`$, follow from those of $`A`$ by taking $`\lambda \nu `$ and $`\lambda _0\nu _0`$. Henceforth, once we have checked the boundary expressions for $`A`$, those for $`C`$ can be checked in a similar way.

Upon substitution of $`\lambda =\lambda _0`$ in the expression for $`A`$ (4.4.1), the first integrand is proportional to $`z_1\lambda ^{}`$ and thus is analytic within the contour $`C_1^\lambda `$. The contribution to the boundary expression therefore needs to come from either $`C_1^\mu `$ or $`C_1^\nu `$. The substitution

$$z_1\lambda _0=\frac{\lambda _0\nu }{\mu \nu }(z_1\mu )\frac{\lambda _0\mu }{\mu \nu }(z_1\nu ),$$
(4.24)
splits the first integral into two complete elliptic integrals

$$\begin{array}{c}\frac{\lambda _0\nu }{\mu \nu }\underset{C_1}{}\frac{(z_1\mu )^{\frac{1}{2}}\mathrm{d}z_1}{(z_1\nu )^{\frac{1}{2}}(z_1\mu _0)^{\frac{1}{2}}(z_1\nu _0)^{\frac{3}{2}}}\hfill \\ \hfill \frac{\lambda _0\mu }{\mu \nu }\underset{C_1}{}\frac{(z_1\nu )^{\frac{1}{2}}\mathrm{d}z_1}{(z_1\mu )^{\frac{1}{2}}(z_1\mu _0)^{\frac{1}{2}}(z_1\nu _0)^{\frac{3}{2}}}.\end{array}$$
(4.25)
Within the contour $`C_1^\mu `$, the integrals can be evaluated in terms of $`E^{}(u)`$ and $`E(u)`$ respectively. When $`\lambda =\lambda _0`$, the second integral in (4.4.1) has a single pole contribution from the contour $`C_2^\lambda `$. Together, $`C_1^\mu C_2^\lambda `$, exactly reproduces the boundary expression $`A(\lambda _0,\mu ,\nu )`$ in (4.2).

When $`\mu =\mu _0`$, both integrands in the expression for $`A`$ have a single pole within the contour $`C_i^\mu `$. However, the combination $`C_1^\mu C_2^\mu `$ does not give the correct boundary expression. We again split both integrals to obtain the required complete elliptic integrals. In the first we substitute

$$z_1\lambda _0=\frac{\lambda _0\nu _0}{\mu _0\nu _0}(z_1\mu _0)\frac{\lambda _0\mu _0}{\mu _0\nu _0}(z_1\nu _0).$$
(4.26)
For the contour $`C_1^\lambda `$, the first integral after the split can be evaluated in terms of $`E^{}(v)`$, but the second integral we leave unchanged. For the integral in the $`z_2`$-plane we substitute

$$z_2\nu _0=\frac{\lambda _0\nu _0}{\lambda _0\mu _0}(z_2\mu _0)\frac{\mu _0\nu _0}{\lambda _0\mu _0}(z_2\lambda _0).$$
(4.27)
We take $`C_2^\nu `$ as contour, and evaluate the first integral after the split in terms of $`E(v)`$. We again leave the second integral unchanged. Except for the contour choice, it is of the same form as the integral we left unchanged in the $`z_1`$-plane.

To obtain the required boundary expression for $`A`$ at $`\mu =\mu _0`$, it turns out that we have to add the contribution of three pairs of contours, $`C_1^\lambda C_2^\mu `$, $`C_1^\mu C_2^\nu `$ and $`C_1^\mu C_2^\mu `$. With the above substitutions (4.26) and (4.27), the first two pairs together provide the required boundary expression, but in addition we have two similar contour integrals

$$\frac{i/8\pi }{(\lambda _0\nu _0)^{\frac{1}{2}}(\lambda \nu )^{\frac{1}{2}}}\underset{C^\tau }{}\frac{(z\lambda )^{\frac{1}{2}}\mathrm{d}z}{(z\nu )^{\frac{1}{2}}(z\lambda _0)^{\frac{1}{2}}(z\nu _0)^{\frac{1}{2}}(z\mu _0)},$$
(4.28)
with contours $`C^\lambda `$ and $`C^\nu `$, respectively. The third pair, $`C_1^\mu C_2^\mu `$, involves the product of two single pole contributions. The resulting polynomial

$$\frac{i/8\pi }{(\lambda _0\nu _0)^{\frac{1}{2}}(\lambda \nu )^{\frac{1}{2}}}\frac{2\pi i(\lambda \mu _0)^{\frac{1}{2}}}{(\mu _0\nu )^{\frac{1}{2}}(\lambda _0\mu _0)^{\frac{1}{2}}(\mu _0\nu _0)^{\frac{1}{2}}},$$
(4.29)
can be written in the same form as (4.28), with contour $`C^\mu `$. As a result, we now have the same integral over all three contours, so that from (4.23), the cumulative result vanishes and we are left with the required boundary expression.

The expression for $`A`$ at $`\nu =\nu _0`$ resembles the one for $`B`$ at the same boundary. This is expected since their boundary expressions in (4.2) are also very similar. The first integral now has a contribution from a double pole in the contour $`C_1^\nu `$. The second integral has no contribution from the contour $`C_2^\nu `$. However, within $`C_2^\mu `$, the second integral can be evaluated in terms of $`E(w)`$. We obtain the correct boundary expression $`A(\lambda ,\mu ,\nu _0)`$ by considering the pair $`C_1^\nu C_2^\mu `$.

#### 4.4.3 Combination of contours

In the previous paragraphs we have constructed a homogeneous solution for $`A`$, $`B`$ and $`C`$, and we have shown that with this solution all nine boundary expressions can be satisfied. For each boundary expression separately, we have determined the required pair of contours and also contours from which there is no contribution. Now we have to find the right combination of all these contours to fit the boundary expressions simultaneously.

We first summarise the required and non-contributing pairs of contours per boundary expression

$`A(\lambda _0,\mu ,\nu )`$ $`:`$ $`C_1^\mu C_2^\lambda \pm C_1^\lambda C_2^\tau ,`$
$`A(\lambda ,\mu _0,\nu )`$ $`:`$ $`+C_1^\mu C_2^\nu +C_1^\lambda C_2^\mu +C_1^\mu C_2^\mu ,`$
$`A(\lambda ,\mu ,\nu _0)`$ $`:`$ $`C_1^\nu C_2^\mu \pm C_1^\tau C_2^\nu ,`$
$`B(\lambda _0,\mu ,\nu )`$ $`:`$ $`C_1^\mu C_2^\lambda \pm C_1^\lambda C_2^\tau ,`$
$`B(\lambda ,\mu _0,\nu )`$ $`:`$ $`\pm C_1^\mu C_2^\tau \pm C_1^\tau C_2^\mu ,`$ (4.30)
$`B(\lambda ,\mu ,\nu _0)`$ $`:`$ $`C_1^\nu C_2^\mu \pm C_1^\tau C_2^\nu ,`$
$`C(\lambda _0,\mu ,\nu )`$ $`:`$ $`C_1^\mu C_2^\lambda \pm C_1^\lambda C_2^\tau ,`$
$`C(\lambda ,\mu _0,\nu )`$ $`:`$ $`+C_1^\mu C_2^\nu +C_1^\lambda C_2^\mu +C_1^\mu C_2^\mu ,`$
$`C(\lambda ,\mu ,\nu _0)`$ $`:`$ $`C_1^\nu C_2^\mu \pm C_1^\tau C_2^\nu ,`$
where $`\tau `$ can be $`\lambda `$, $`\mu `$ or $`\nu `$. At each boundary separately, $`\lambda =\lambda _0`$, $`\mu =\mu _0`$ and $`\nu =\nu _0`$, the allowed combination of contours matches between $`A`$, $`B`$ and $`C`$. This leaves the question how to relate the combination of contours at the different boundaries relative to each other.

From (4.23), we know that in both the complex $`z_1`$-plane and $`z_2`$-plane, the cumulative contribution of the three simple contours cancels. As a consequence, each of the following three combinations of integration contours

$$C_1^\mu C_2^\mu =C_1^\mu (C_2^\lambda +C_2^\nu )=(C_1^\lambda +C_1^\nu )C_2^\mu ,$$
(4.31)
will give the same result. Similarly, we can add to each combination the pairs $`C_1^\lambda C_2^\mu `$ and $`C_1^\mu C_2^\nu `$, to obtain

$$C_1^\mu C_2^\nu +C_1^\lambda C_2^\mu +C_1^\mu C_2^\mu =C_1^\lambda C_2^\mu C_1^\mu C_2^\lambda =C_1^\mu C_2^\nu C_1^\nu C_2^\mu .$$
(4.32)
The first combination of contour pairs matches the allowed range for $`\mu =\mu _0`$ in (4.4.3) and the second and third match the boundary expressions $`\lambda =\lambda _0`$ and $`\nu =\nu _0`$. This completes the proof that the expressions (4.4.1)–(4.4.1) for $`A`$, $`B`$ and $`C`$ solve the homogeneous equations (4.2) and satisfy the nine boundary expressions (4.2) simultaneously when the integration contour is any of the three combinations (4.32). We shall see below that the first of these combinations is preferred in numerical evaluations.

### 4.5 Evaluation of the homogeneous solutions

We write the complex contour integrals in the homogeneous solutions $`A`$, $`B`$ and $`C`$ (4.4.14.4.1) as real integrals. The resulting complete hyperelliptic integrals are expressed as single quadratures, which can be evaluated numerically in a straightforward way. We also express the complete elliptic integrals in the two-dimensional homogeneous solutions $`F`$, $`G`$, $`H`$ and $`I`$ (4.2) in this way to facilitate their numerical evaluation.

#### 4.5.1 From complex to real integrals

To transform the complex contour integrals in (4.4.1)–(4.4.1) in real integrals we wrap the contours $`C^\lambda `$, $`C^\mu `$ and $`C^\nu `$ around the corresponding pair of branch points (Fig. 6). The integrands consist of terms $`z\tau `$ and $`z\tau _0`$, all with powers larger than $`1`$, except $`z_1\nu _0`$ and $`z_2\lambda _0`$, both of which occur to the power $`\frac{3}{2}`$. This means that for all simple contours $`C_i^\tau `$ $`(\tau =\lambda ,\mu ,\nu ;i=1,2)`$, except for $`C_1^\nu `$ and $`C_2^\lambda `$, the contribution from the arcs around the branch points vanishes. In the latter case, we are left with the parts parallel to the real axis, so that we can rewrite the complex integrals as real integrals with the branch points as integration limits. The only combination of contours of the three given in (4.32) that does not involve both $`C_1^\nu `$ and $`C_2^\lambda `$, is

$$SC_1^\mu C_2^\nu +C_1^\lambda C_2^\mu +C_1^\mu C_2^\mu .$$
(4.33)
We have to be careful with the changes in phase when wrapping each of the simple contours around the branch points. One can verify that the phase changes per contour are the same for all three the homogeneous solutions $`A`$, $`B`$ and $`C`$, and also that the contribution from the parts parallel to the real axis is equivalent. This gives a factor 2 per contour and thus a factor 4 for the combination of contour pairs in $`S`$. In this way, we can transform the double complex contour integration into the following combination of real integrals

$$\underset{S}{}dz_1dz_2=4(\underset{\lambda }{\overset{\lambda _0}{}}dt_1\underset{\mu }{\overset{\mu _0}{}}dt_2+\underset{\mu }{\overset{\mu _0}{}}dt_1\underset{\nu }{\overset{\nu _0}{}}dt_2\underset{\mu }{\overset{\mu _0}{}}dt_1\underset{\mu }{\overset{\mu _0}{}}dt_2),$$
(4.34)
with $`t_i`$ the real part of $`z_i`$.

We apply this transformation to (4.4.1)–(4.4.1), and we absorb the factor of 4 left in the denominators into the integrals, so that we can write

$`A(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)`$ $`=`$ $`{\displaystyle \frac{(\mu _0\nu _0)(\mu _0\lambda _0)\mathrm{\Lambda }}{\pi ^2(\lambda \mu )(\lambda \nu )}}(A_1A_2+A_3A_4A_2A_3),`$
$`B(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)`$ $`=`$ $`{\displaystyle \frac{(\mu _0\nu _0)(\mu _0\lambda _0)\mathrm{\Lambda }}{\pi ^2(\mu \nu )(\mu \lambda )}}(B_1B_2+B_3B_4B_2B_3),`$
$`C(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)`$ $`=`$ $`{\displaystyle \frac{(\mu _0\nu _0)(\mu _0\lambda _0)\mathrm{\Lambda }}{\pi ^2(\nu \lambda )(\nu \mu )}}(C_1C_2+C_3C_4C_2C_3),`$
where $`A_i`$, $`B_i`$ and $`C_i`$ ($`i=1,2,3,4`$) are complete hyperelliptic integrals, for which we give expressions below, and

$$\mathrm{\Lambda }^2=\frac{(\lambda \mu )(\lambda \nu )(\mu \nu )}{(\lambda _0\mu _0)(\lambda _0\nu _0)(\mu _0\nu _0)}.$$
(4.36)
The second set of arguments added to $`A`$, $`B`$ and $`C`$ make explicit the position $`(\lambda _0,\mu _0,\nu _0)`$ of the source point which is causing the stresses at the field point $`(\lambda ,\mu ,\nu )`$.

#### 4.5.2 The complete hyperelliptic integrals

With the transformation described in the previous section the expression for, e.g., the complete hyperelliptic integral $`A_2`$ is of the form

$$A_2=\frac{1}{2}\underset{\mu }{\overset{\mu _0}{}}\frac{\mathrm{d}t}{\lambda _0t}\sqrt{\frac{(\lambda t)(t\nu _0)}{(\mu _0t)(t\mu )(\lambda _0t)(t\nu )}}.$$
(4.37)
The integrand has two singularities, one at the lower integration limit $`t=\mu `$ and one at the upper integration limit $`t=\mu _0`$. The substitution $`t=\mu +(\mu _0\mu )\mathrm{cos}^2\theta `$ removes both singularities, since $`\mathrm{d}t/\sqrt{(\mu _0t)(t\mu )}=2(\mu _0\mu )\mathrm{d}\theta `$.

All complete hyperelliptic integrals $`A_i`$, $`B_i`$ and $`C_i`$ ($`i=1,2,3,4`$) in (4.5.1) are of the form (4.37) and have at most two singularities at either of the integration limits. Hence, we can apply a similar substitution to remove the singularities. This results in the following expressions

$`A_1`$ $`=(\lambda _0\lambda )^2{\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{\mathrm{sin}^2\theta \mathrm{cos}^2\theta \mathrm{d}\theta }{x_3\mathrm{\Delta }_x}},`$ $`A_2`$ $`={\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{y_1y_4\mathrm{d}\theta }{y_3\mathrm{\Delta }_y}},`$
$`A_4`$ $`=(\nu _0\nu ){\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{z_2\mathrm{sin}^2\theta \mathrm{d}\theta }{z_1\mathrm{\Delta }_z}},`$ $`A_3`$ $`={\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{y_3y_4\mathrm{d}\theta }{y_1\mathrm{\Delta }_y}},`$
$`B_1`$ $`=(\lambda _0\lambda ){\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{x_2\mathrm{sin}^2\theta \mathrm{d}\theta }{x_3\mathrm{\Delta }_x}},`$ $`B_2`$ $`=(\mu _0\mu ){\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{y_1\mathrm{cos}^2\theta \mathrm{d}\theta }{y_3\mathrm{\Delta }_y}},`$
$`B_4`$ $`=(\nu _0\nu ){\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{z_4\mathrm{sin}^2\theta \mathrm{d}\theta }{z_1\mathrm{\Delta }_z}},`$ $`B_3`$ $`=(\mu _0\mu ){\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{y_3\mathrm{cos}^2\theta \mathrm{d}\theta }{y_1\mathrm{\Delta }_y}},`$
$`C_1`$ $`=(\lambda _0\lambda ){\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{x_4\mathrm{sin}^2\theta \mathrm{d}\theta }{x_3\mathrm{\Delta }_x}},`$ $`C_2`$ $`={\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{y_1y_2\mathrm{d}\theta }{y_3\mathrm{\Delta }_y}},`$
$`C_4`$ $`=(\nu _0\nu )^2{\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{\mathrm{sin}^2\theta \mathrm{cos}^2\theta \mathrm{d}\theta }{z_1\mathrm{\Delta }_z}},`$ $`C_3`$ $`={\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{y_2y_3\mathrm{d}\theta }{y_1\mathrm{\Delta }_y}},`$
where we have defined

$$\mathrm{\Delta }_x^2=x_1x_2x_3x_4,\mathrm{\Delta }_y^2=y_1y_2y_3y_4,\mathrm{\Delta }_z^2=z_1z_2z_3z_4,$$
(4.39)
and the factors $`x_i`$, $`y_i`$ and $`z_i`$ $`(i=1,2,3,4)`$ are given by

$`x_1`$ $`=(\lambda \mu _0)+(\lambda _0\lambda )\mathrm{cos}^2\theta ,`$ $`x_2`$ $`=(\lambda \mu )+(\lambda _0\lambda )\mathrm{cos}^2\theta ,`$
$`x_3`$ $`=(\lambda \nu _0)+(\lambda _0\lambda )\mathrm{cos}^2\theta ,`$ $`x_4`$ $`=(\lambda \nu )+(\lambda _0\lambda )\mathrm{cos}^2\theta ,`$
$`y_1`$ $`=(\mu \nu _0)+(\mu _0\mu )\mathrm{cos}^2\theta ,`$ $`y_2`$ $`=(\mu \nu )+(\mu _0\mu )\mathrm{cos}^2\theta ,`$
$`y_3`$ $`=(\mu \lambda _0)+(\mu _0\mu )\mathrm{cos}^2\theta ,`$ $`y_4`$ $`=(\mu \lambda )+(\mu _0\mu )\mathrm{cos}^2\theta ,`$
$`z_1`$ $`=(\nu \lambda _0)+(\nu _0\nu )\mathrm{cos}^2\theta ,`$ $`z_2`$ $`=(\nu \lambda )+(\nu _0\nu )\mathrm{cos}^2\theta ,`$
$`z_3`$ $`=(\nu \mu _0)+(\nu _0\nu )\mathrm{cos}^2\theta ,`$ $`z_4`$ $`=(\nu \mu )+(\nu _0\nu )\mathrm{cos}^2\theta .`$ (4.40)
For each $`i`$ these factors follow from each other by cyclic permutation of $`\lambda \mu \nu \lambda `$ and at the same time $`\lambda _0\mu _0\nu _0\lambda _0`$. Half of the factors – all $`x_i`$, $`y_1`$ and $`y_2`$ – are always positive, whereas the other factors are always negative. The latter implies that one has to be careful with the signs of the factors under the square root when evaluating the single quadratures numerically.

#### 4.5.3 The complete elliptic integrals

The two-dimensional homogeneous solutions $`F`$, $`G`$, $`H`$ and $`I`$ are given in (4.2) in terms of the Legendre complete elliptic integrals $`E(m)`$ and $`E^{}(m)=[E(m)K(m)]/2m`$. Numerical routines for $`E(m)`$ and $`K(m)`$ (e.g., Press et al. 1999) generally require the argument to be $`0m<1`$. In the allowed range of the confocal ellipsoidal coordinates, the arguments $`u`$ (4.10) and $`w`$ (3.16) become larger than unity. In these cases we can use transformations to express $`E(m)`$ and $`K(m)`$ in terms of $`E(1/m)`$ and $`K(1/m)`$ (e.g., Byrd & Friedman 1971).

We prefer, however, to write the complete elliptic integrals as single quadratures similar to the above expressions for the hyperelliptic integrals. These quadratures can easily be evaluated numerically and apply to the full range of the confocal ellipsoidal coordinates. The resulting expressions for the two-dimensional homogeneous solutions are

$`F(\lambda ,\mu ;\lambda _0,\mu _0)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\sqrt{{\displaystyle \frac{\lambda \mu }{\lambda _0\mu _0}}}{\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{x_1d\theta }{x_2\sqrt{x_1x_2}}},`$
$`G(\lambda ,\mu ;\lambda _0,\mu _0)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\sqrt{{\displaystyle \frac{\lambda \mu }{\lambda _0\mu _0}}}(\mu _0\mu ){\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{\mathrm{sin}^2\theta d\theta }{y_4\sqrt{y_3y_4}}},`$
$`H(\mu ,\nu ;\mu _0,\nu _0)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\sqrt{{\displaystyle \frac{\mu \nu }{\mu _0\nu _0}}}(\mu _0\mu ){\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{\mathrm{sin}^2\theta d\theta }{y_2\sqrt{y_1y_2}}},`$
$`I(\mu ,\nu ;\mu _0,\nu _0)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\sqrt{{\displaystyle \frac{\mu \nu }{\mu _0\nu _0}}}{\displaystyle \underset{0}{\overset{\pi /2}{}}}{\displaystyle \frac{z_3d\theta }{z_4\sqrt{z_3z_4}}}.`$ (4.41)
Again we have added two arguments to make the position of the unit source explicitly. We note that the homogeneous solutions $`A(\lambda ,\mu ;\lambda _0,\mu _0)`$ and $`B(\lambda ,\mu ;\lambda _0,\mu _0)`$ for the disc case (3.41) are equivalent to $`F`$ and $`G`$ respectively.

### 4.6 General triaxial solution

We now construct the solution of the Jeans equations for triaxial Stäckel models (4.2), by superposition of singular solutions, which involve the homogeneous solution derived in the above. We match the solution to the boundary conditions at $`\mu =\alpha `$ and $`\nu =\beta `$, and check for convergence of the solution when $`\lambda \mathrm{}`$. Next, we consider alternative boundary conditions and present the triaxial solution for a general finite region. We also show that the general solution yields the correct result in the case of thin tube orbits and the triaxial Abel models of Dejonghe & Laurent (1991). Finally, we describe a numerical test of the triaxial solution to a polytrope model.

#### 4.6.1 Superposition of singular solutions

Substitution of the functions $`A`$, $`B`$, $`C`$ (4.5.1) and the functions $`F`$, $`G`$, $`H`$, $`I`$ (4.5.3) in expression (4.2), provides the three singular solutions of the system of simplified Jeans equations, with the right-hand side given by (4.1). We denote these by $`S_2^{\tau \tau }`$ $`(\tau =\lambda ,\mu ,\nu )`$. The singular solutions of the two similar simplified systems, with the triplet of delta functions at the right-hand side of the first and third equation, $`S_1^{\tau \tau }`$ and $`S_3^{\tau \tau }`$ then follow from $`S_2^{\tau \tau }`$ by cyclic permutation. This gives

$`S_1^{\lambda \lambda }`$ $`=`$ $`B(\nu ,\lambda ,\mu ;\nu _0,\lambda _0,\mu _0)+G(\nu ,\lambda ;\nu _0,\lambda _0)\delta (\mu _0\mu )`$
$`+H(\lambda ,\mu ;\lambda _0,\mu _0)\delta (\nu _0\nu )\delta (\mu _0\mu )\delta (\nu _0\nu ),`$
$`S_1^{\mu \mu }`$ $`=`$ $`C(\nu ,\lambda ,\mu ;\nu _0,\lambda _0,\mu _0)+I(\lambda ,\mu ;\lambda _0,\mu _0)\delta (\nu _0\nu )`$
$`S_1^{\nu \nu }`$ $`=`$ $`A(\nu ,\lambda ,\mu ;\nu _0,\lambda _0,\mu _0)+F(\nu ,\lambda ;\nu _0,\lambda _0)\delta (\mu _0\mu ),`$ (4.42a)
$`S_2^{\lambda \lambda }`$ $`=`$ $`A(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)+F(\lambda ,\mu ;\lambda _0,\mu _0)\delta (\nu _0\nu ),`$
$`S_2^{\mu \mu }`$ $`=`$ $`B(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)+G(\lambda ,\mu ;\lambda _0,\mu _0)\delta (\nu _0\nu )`$
$`+H(\mu ,\nu ;\mu _0,\nu _0)\delta (\lambda _0\lambda )\delta (\nu _0\nu )\delta (\lambda _0\lambda ),`$
$`S_2^{\nu \nu }`$ $`=`$ $`C(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)+I(\mu ,\nu ;\mu _0,\nu _0)\delta (\lambda _0\lambda )`$ (4.42b)
$`S_3^{\lambda \lambda }`$ $`=`$ $`C(\mu ,\nu ,\lambda ;\mu _0,\nu _0,\lambda _0)+I(\nu ,\lambda ;\nu _0,\lambda _0)\delta (\mu _0\mu ),`$
$`S_3^{\mu \mu }`$ $`=`$ $`A(\mu ,\nu ,\lambda ;\mu _0,\nu _0,\lambda _0)+F(\mu ,\nu ;\mu _0,\nu _0)\delta (\lambda _0\lambda )`$
$`S_3^{\nu \nu }`$ $`=`$ $`B(\mu ,\nu ,\lambda ;\mu _0,\nu _0,\lambda _0)+G(\mu ,\nu ;\mu _0,\nu _0)\delta (\lambda _0\lambda )`$ (4.42c)
$`+H(\nu ,\lambda ;\nu _0,\lambda _0)\delta (\mu _0\mu )\delta (\lambda _0\lambda )\delta (\mu _0\mu ).`$
These singular solutions describe the contribution of a source point in $`(\lambda _0,\mu _0,\nu _0)`$ to $`(\lambda ,\mu ,\nu )`$. To find the solution of the full equations (4.2), we multiply the singular solutions (4.42), (4.42) and (4.42) by $`g_1(\lambda _0,\mu _0,\nu _0)`$, $`g_2(\lambda _0,\mu _0,\nu _0)`$ and $`g_3(\lambda _0,\mu _0,\nu _0)`$, respectively, so that the contribution from the source point naturally depends on the local density and potential (cf. eq. \[4.3\]). Then, for each coordinate $`\tau =\lambda ,\mu ,\nu `$, we add the three weighted singular solutions, and integrate over the volume $`\mathrm{\Omega }`$, defined as

$$\mathrm{\Omega }=\{(\lambda _0,\mu _0,\mu _0):\lambda \lambda _0<\mathrm{},\mu \mu _0\alpha ,\nu \nu _0\beta \},$$
(4.43)
which is the three-dimensional extension of the integration domain $`D`$ in Fig. 4. The resulting solution solves the inhomogeneous Jeans equations (4.2), but does not give the correct values at the boundaries $`\mu =\alpha `$ and $`\nu =\beta `$. They are found by multiplying the singular solutions (4.42) evaluated at $`\mu _0=\alpha `$, and, similarly, the singular solutions (4.42) evaluated at $`\nu _0=\beta `$, by $`S_{\mu \mu }(\lambda _0,\alpha ,\nu _0)`$ and $`S_{\nu \nu }(\lambda _0,\mu _0,\beta )`$, respectively, and integrating in $`\mathrm{\Omega }`$ over the coordinates that are not fixed. One can verify that this procedure represents the boundary values correctly. The final result for the general solution of the Jeans equations (4.2) for triaxial Stäckel models is

$`S_{\tau \tau }(\lambda ,`$ $`\mu `$ $`,\nu )={\displaystyle }_\lambda ^{\mathrm{}}\mathrm{d}\lambda _0{\displaystyle }_\mu ^\alpha \mathrm{d}\mu _0{\displaystyle }_\nu ^\beta \mathrm{d}\nu _0{\displaystyle }_{i=1}^3g_i(\lambda _0,\mu _0,\nu _0)S_i^{\tau \tau }(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)`$ (4.44)
$``$ $`{\displaystyle \underset{\nu }{\overset{\beta }{}}}d\nu _0{\displaystyle \underset{\lambda }{\overset{\mathrm{}}{}}}d\lambda _0S_{\mu \mu }(\lambda _0,\alpha ,\nu _0)S_2^{\tau \tau }(\lambda ,\mu ,\nu ;\lambda _0,\alpha ,\nu _0)`$
$``$ $`{\displaystyle \underset{\lambda }{\overset{\mathrm{}}{}}}d\lambda _0{\displaystyle \underset{\mu }{\overset{\alpha }{}}}d\mu _0S_{\nu \nu }(\lambda _0,\mu _0,\beta )S_3^{\tau \tau }(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\beta ),`$
where $`\tau =(\lambda ,\mu ,\nu )`$. This gives the stresses everywhere, once we have specified $`S_{\mu \mu }(\lambda ,\alpha ,\nu )`$ and $`S_{\nu \nu }(\lambda ,\mu ,\beta )`$. At both boundaries $`\mu =\alpha `$ and $`\nu =\beta `$, the three stress components are related by a set of two Jeans equations, i.e., (4.2) evaluated at $`\mu =\alpha `$ and $`\nu =\beta `$ respectively. From §3, we know that the solution of these two-dimensional systems both will involve a (boundary) function of one variable. We need this latter freedom to satisfy the continuity conditions (2.17). This means it is sufficient to specify any of the three stress components at $`\mu =\alpha `$ and $`\nu =\beta `$.

#### 4.6.2 Convergence of the general triaxial solution

As in §§3.1.4, 3.2.7 and 3.4 we suppose $`G(\tau )=𝒪(\tau ^\delta )`$ when $`\tau \mathrm{}`$, with $`\delta `$ in the range $`[\frac{1}{2},0)`$. This implies that the potential $`V_S`$ (2.3) is also $`𝒪(\tau ^\delta )`$. We assume that the density $`\rho `$, which does not need to be the density $`\rho _S`$ which generates $`V_S`$, is of the form $`N(\mu ,\nu )\lambda ^{s/2}`$ when $`\lambda \mathrm{}`$. In the special case where $`\rho =\rho _S`$, we have $`s4`$ except possibly along the $`z`$-axis. When $`s=4`$ the models remain flattened out to the largest radii, but when $`s<4`$ the function $`N(\mu ,\nu )1`$ in the limit $`\lambda \mathrm{}`$ (de Zeeuw et al. 1986).

From the definition (4.3), we find that $`g_1(\lambda _0,\mu _0,\nu _0)=𝒪(\lambda _0^{\delta s/2})`$ as $`\lambda _0\mathrm{}`$, while $`g_2(\lambda _0,\mu _0,\nu _0)`$ and $`g_3(\lambda _0,\mu _0,\nu _0)`$ are larger and both $`𝒪(\lambda _0^{s/2})`$. To investigate the behaviour of the singular solutions (4.42) at large distance, we have to carefully analyse the complete hyperelliptic (4.38) and elliptic (4.5.3) integrals as $`\lambda _0\mathrm{}`$. This is simplified by writing them as Carlson’s $`R`$-functions (Carlson 1977). We finally find for the singular solutions that $`S_1^{\tau \tau }=𝒪(1)`$ when $`\lambda _0\mathrm{}`$, whereas $`S_2^{\tau \tau }`$ and $`S_3^{\tau \tau }`$ are smaller and $`𝒪(\lambda _0^1)`$, with $`\tau =\lambda ,\mu ,\nu `$. This shows that for the volume integral in the triaxial solution (4.44) to converge, we must have $`\delta s/2+1<0`$. This is equivalent to the requirement $`s>2\delta +2`$ we obtained in §3.4 for the limiting cases of prolate and oblate potentials and for the large radii limit with scale-free DF. From the convergence of the remaining two double integrals in (4.44), we find that the boundary stresses $`S_{\mu \mu }(\lambda ,\alpha ,\nu )`$ and $`S_{\nu \nu }(\lambda ,\mu ,\beta )`$ cannot exceed $`𝒪(1)`$ when $`\lambda \mathrm{}`$.

The latter is in agreement with the large $`\lambda `$ behaviour of $`S_{\tau \tau }(\lambda ,\mu ,\nu )`$ that follows from the volume integral. The singular solutions $`S_i^{\lambda \lambda }=𝒪(1)`$ ($`i=1,2,3`$) when $`\lambda \mathrm{}`$, larger than $`S_i^{\mu \mu }`$ and $`S_i^{\nu \nu }`$, which are all $`𝒪(\lambda ^1)`$. Evaluating the volume integral at large distance then gives $`S_{\tau \tau }(\lambda ,\mu ,\nu )=𝒪(\lambda ^{\delta s/2+1})`$, i.e., not exceeding $`𝒪(1)`$ if the requirement $`s>2\delta +2`$ is satisfied. We obtain the same behaviour and requirement from the energy equation (2.10).

We conclude that for the general triaxial case, as well as for the limiting cases with a three-dimensional shape, the stress components $`T_{\tau \tau }(\lambda ,\mu ,\nu )`$ are $`𝒪(\lambda ^{\delta s/2})`$ at large distance, with the requirement that $`s>2\delta +2`$ for $`\frac{1}{2}\delta <0`$. We obtained the same result for the stresses in the disc case, except that then $`s>2\delta +1`$. Both the three-dimensional and two-dimensional requirements are met for many density distributions $`\rho `$ and potentials $`V_S`$ of interest. They do not break down until the isothermal limit $`\delta 0`$, with $`s=1`$ (disc) and $`s=2`$ (three-dimensional) is reached.

#### 4.6.3 Alternative boundary conditions

Our solution for the stress components at each point $`(\lambda ,\mu ,\nu )`$ in a triaxial model with a Stäckel potential consists of the weighted contribution of all sources outwards of this point. Accordingly, we have integrated with respect to $`\lambda _0`$, $`\mu _0`$ and $`\nu _0`$, with lower limits the coordinates of the chosen point and upper limits $`\mathrm{}`$, $`\alpha `$ and $`\beta `$, respectively. To obtain the correct expressions at the outer boundaries, the stresses must vanish when $`\lambda \mathrm{}`$ and they have to be specified at $`\mu =\alpha `$ and $`\nu =\beta `$.

The integration limits $`\lambda `$, $`\mu `$ and $`\nu `$ are fixed, but for the other three limits we can, in principle, equally well choose $`\alpha `$, $`\beta `$ and $`\gamma `$ respectively. The latter choices also imply the specification of the stress components at these boundaries instead. Each of the eight possible combinations of these limits corresponds to one of the octants into which the physical region $`\gamma \nu _0\beta \mu _0\alpha \lambda _0<\mathrm{}`$ is split by the lines through the point $`(\lambda ,\mu ,\nu )`$. By arguments similar to those given in §3.3, one may show that in all octants the expressions (4.5.1) for $`A`$, $`B`$, $`C`$, and (4.2) for $`F`$, $`G`$, $`H`$, $`I`$ are equivalent. Hence, again the only differences in the singular solutions are due to possible changes in the sign of the step-functions, but the changes in the integration limits cancel the sign differences between the corresponding singular solutions. However, as in §3.3 for the two-dimensional case, it is not difficult to show that while switching the boundary conditions $`\mu `$ and $`\nu `$ is indeed straightforward, the switch from $`\lambda \mathrm{}`$ to $`\lambda =\alpha `$ again leads to solutions which generally have the incorrect radial fall-off, and hence are non-physical.

#### 4.6.4 Triaxial solution for a general finite region

If we denote non-fixed integration limits by $`\lambda _e`$, $`\mu _e`$ and $`\nu _e`$ respectively, we can write the triaxial solution for a general finite region as

$`S_{\tau \tau }(\lambda ,`$ $`\mu `$ $`,\nu )={\displaystyle }_\lambda ^{\lambda _e}\mathrm{d}\lambda _0{\displaystyle }_\mu ^{\mu _e}\mathrm{d}\mu _0{\displaystyle }_\nu ^{\nu _e}\mathrm{d}\nu _0{\displaystyle }_{i=1}^3g_i(\lambda _0,\mu _0,\nu _0)S_i^{\tau \tau }(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)`$ (4.45)
$``$ $`{\displaystyle \underset{\mu }{\overset{\mu _e}{}}}d\mu _0{\displaystyle \underset{\nu }{\overset{\nu _e}{}}}d\nu _0S_{\lambda \lambda }(\lambda _e,\mu _0,\nu _0)S_1^{\tau \tau }(\lambda ,\mu ,\nu ;\lambda _e,\mu _0,\nu _0)`$
$``$ $`{\displaystyle \underset{\nu }{\overset{\nu _e}{}}}d\nu _0{\displaystyle \underset{\lambda }{\overset{\lambda _e}{}}}d\lambda _0S_{\mu \mu }(\lambda _0,\mu _e,\nu _0)S_2^{\tau \tau }(\lambda ,\mu ,\nu ;\lambda _0,\mu _e,\nu _0)`$
$``$ $`{\displaystyle \underset{\lambda }{\overset{\lambda _e}{}}}d\lambda _0{\displaystyle \underset{\mu }{\overset{\mu _e}{}}}d\mu _0S_{\nu \nu }(\lambda _0,\mu _0,\nu _e)S_3^{\tau \tau }(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _e),`$
with, as usual, $`\tau =\lambda ,\mu ,\nu `$. The weight functions $`g_i`$ ($`i=1,2,3`$) are defined in (4.3) and the singular solutions $`S_i^{\tau \tau }`$ are given by (4.42). The non-fixed integration limits are chosen in the corresponding physical ranges, i.e., $`\lambda _e[\alpha ,\mathrm{}]`$, $`\mu _e[\beta ,\alpha ]`$ and $`\nu _e[\gamma ,\beta ]`$, but $`\lambda _e\alpha `$ (see §4.6.3). The solution requires the specification of the stress components on the boundary surfaces $`\lambda =\lambda _e`$, $`\mu =\mu _e`$ and $`\nu =\nu _e`$. On each of these surfaces the three stress components are related by two of the three Jeans equations (4.2) and the continuity conditions (2.17). Hence, once one of the stress components is prescribed on three boundary surfaces, the solution (4.44) yields all three stresses everywhere in the triaxial Stäckel galaxy. The stresses on the remaining three boundary surfaces then follow as the limits of the latter solution.

#### 4.6.5 Physical solutions

Statler (1987) and HZ92 showed that many different DFs are consistent with a triaxial density $`\rho `$ in the potential $`V_S`$. Specifically, the boundary plane $`\nu =\beta `$, i.e., the area outside the focal hyperbola in the ($`x,z`$)-plane (Fig. 2), is only reached by inner (I) and outer (O) long-axis tube orbits. A split between the contribution of both orbit families to the density in this plane has to be chosen, upon which the DF for both the I and O orbits is fixed in case only thin tubes are populated, but many other possibilities exist when the full set of I- and O-orbits is included. For each of these DFs, the density provided by the I- and O-tubes can then in principle be found throughout configuration space. In the area outside the focal ellipse in the ($`y,z`$)-plane ($`\mu =\alpha `$), only the O-tubes and S-tubes contribute to the density. Subtracting the known density of the O-orbits leaves the density to be provided by the S-tubes in this plane, from which their DF can be determined. This is again unique when only thin orbits are used, but is non-unique otherwise. The density that remains after subtracting the I-, O-, and S-tube densities from $`\rho `$ must be provided by the box (B) orbits. Their DF is now fixed, and can be found by solving a system of linear equations, starting from the outside ($`\lambda \mathrm{}`$).

The total DF is the sum of the DFs of the four orbit families, and is hence highly non-unique. All these DFs give rise to a range of stresses $`T_{\lambda \lambda },T_{\mu \mu },T_{\nu \nu }`$, and our solution of the Jeans equations must be sufficiently general to contain them as a subset. This is indeed the case, as we are allowed to choose the stresses on the special surfaces $`\nu =\beta `$ and $`\mu =\alpha `$. However, not all choices will correspond to physical DFs. The requirement $`T_{\tau \tau }0`$ is necessary but not sufficient for the associated DF to be non-negative everywhere.

#### 4.6.6 The general solution for thin tube orbits

For each of the three tube families in case of infinitesimally thin orbits one of the three stress components vanishes everywhere (see §2.5.6). We are left with two non-zero stress components related to the density and potential by three reduced Jeans equations (4.2). We thus have subsidiary conditions on the three right hand side terms $`g_1`$, $`g_2`$ and $`g_3`$.

HZ92 solved for the two non-trivial stresses and showed that they can be found by single quadratures (with integrands involving no worse than complete elliptic integrals), once the corresponding stress had been chosen at $`\nu =\beta `$ (for I- and O-tubes) or at $`\mu =\alpha `$ (for S-tubes).

By analogy with the reasoning for the thin tube orbits in the disc case (§3.4.4), we can show that for each of the three tube families in the case of thin orbits the general triaxial solution (4.45) gives the stress components correctly. Consider, e.g., the thin I-tubes, for which $`S_{\mu \mu }0`$. Apply the latter to (4.45), substitute for $`g_1`$, $`g_2`$ and $`g_3`$ the subsidiary conditions that follow from the reduced Jeans equations (4.2) and substitute for the singular solutions the expressions (4.42). After several partial integrations, we use that the homogeneous solutions $`A`$, $`B`$ and $`C`$ solve a homogeneous system similar to (4.2), but now with respect to the source point coordinates $`(\lambda _0,\mu _0,\nu _0)`$

$$\frac{B(\nu ,\lambda ,\mu ;\nu _0,\lambda _0,\mu _0)}{\lambda _0}=\frac{A(\lambda ,\mu ,\nu ;\lambda _0,\mu _0,\nu _0)}{2(\lambda _0\mu _0)}+\frac{C(\mu ,\nu ,\lambda ;\mu _0,\nu _0,\lambda _0)}{2(\lambda _0\nu _0)},$$
(4.46)
where other relations follow by cyclic permutation of $`\lambda \mu \nu \lambda `$ and $`\lambda _0\mu _0\nu _0\lambda _0`$. And similar for the two-dimensional homogeneous solutions $`F`$, $`G`$, $`H`$ and $`I`$ the relations follow from

$`{\displaystyle \frac{G(\mu ,\lambda ;\mu _0,\lambda _0)}{\lambda _0}}`$ $`=`$ $`{\displaystyle \frac{F(\lambda ,\mu ;\lambda _0,\mu _0)}{2(\lambda _0\mu _0)}},`$
$`{\displaystyle \frac{H(\mu ,\nu ;\mu _0,\nu _0)}{\mu _0}}`$ $`=`$ $`{\displaystyle \frac{I(\nu ,\mu ;\nu _0,\mu _0)}{2(\mu _0\nu _0)}}.`$
It indeed turns out that for $`S_{\mu \mu }(\lambda ,\mu ,\nu )`$ all terms cancel on the right hand side of (4.45). The terms that are left in the case of $`S_{\lambda \lambda }`$ and $`S_{\nu \nu }`$ are just eq. (4.2a) integrated with respect to $`\lambda `$ and eq. (4.2c) integrated with respect to $`\nu `$, respectively, and using that $`S_{\mu \mu }0`$. A similar analysis as above shows that also for thin O- and S-tubes – for which $`S_{\lambda \lambda }0`$ in both cases – the general triaxial solution yields the correct result.

#### 4.6.7 Triaxial Abel models

For a galaxy with a triaxial potential of Stäckel form, the DF is a function of the three exact isolating integrals of motion, $`f(𝐱,𝐯)=f(E,I_2,I_3)`$ (see also §2.2). The expressions for $`E`$, $`I_2`$ and $`I_3`$ in terms of the phase-space coordinates ($`𝐱,𝐯`$) can be found in e.g. Z85. We can thus write the velocity moments of the DF as a triple integral over $`E`$, $`I_2`$ and $`I_3`$. Assuming that the DF is function of only one variable

$$SE+wI_2+uI_3,$$
(4.48)
with $`w`$ and $`u`$ constants, Dejonghe & Laurent (1991) show that the triple integration simplifies to a one-dimensional Abel integration over $`S`$. Even though a DF of this form can only describe a self-consistent model in the spherical case (ellipsoidal hypothesis, see, e.g., Eddington 1915), the Jeans equations do not require self-consistency.

The special Abel form results in a simple analytical relation between the three stress components (Dejonghe & Laurent 1991, their eq. \[5.6\])

$$T_{\mu \mu }=T_{\lambda \lambda }a_{\mu \nu }/a_{\lambda \nu },T_{\nu \nu }=T_{\lambda \lambda }a_{\mu \nu }/a_{\mu \lambda },$$
(4.49)
with

$$a_{\sigma \tau }=(\gamma \alpha )+(\sigma +\alpha )(\tau +\alpha )w(\sigma +\gamma )(\tau +\gamma )u,$$
(4.50)
and $`\sigma ,\tau =\lambda ,\mu ,\nu `$. With these relations we find that

$$\frac{T_{\lambda \lambda }T_{\mu \mu }}{\lambda \mu }=\frac{T_{\lambda \lambda }}{a_{\lambda \nu }}\frac{a_{\lambda \nu }}{\lambda },\frac{T_{\lambda \lambda }T_{\nu \nu }}{\lambda \nu }=\frac{T_{\lambda \lambda }}{a_{\lambda \mu }}\frac{a_{\lambda \mu }}{\lambda }.$$
(4.51)
The first Jeans equation (2.16a) now becomes a first-order partial differential equation for $`T_{\lambda \lambda }`$. This equation can be solved in a straightforward way and provides an elegant and simple expression for the radial stress component

$$\begin{array}{c}T_{\lambda \lambda }(\lambda ,\mu ,\nu )=\sqrt{\frac{a_{\lambda _e\mu }a_{\lambda _e\nu }}{a_{\lambda \mu }a_{\lambda \nu }}}T_{\lambda \lambda }(\lambda _e,\mu ,\nu )\hfill \\ \hfill +\underset{\lambda }{\overset{\lambda _e}{}}d\lambda _0\left[\sqrt{\frac{a_{\lambda _0\mu }a_{\lambda _0\nu }}{a_{\lambda \mu }a_{\lambda \nu }}}\rho \frac{V_S}{\lambda _0}\right].\end{array}$$
(4.52)
The expressions for $`T_{\mu \mu }`$ and $`T_{\nu \nu }`$ follow by application of the ratios (4.49). If we let the boundary value $`\lambda _e\mathrm{}`$, the first term on the right-hand side of (4.52) vanishes.

The density $`\rho `$, which does not need to be the density $`\rho _S`$ which generates $`V_S`$, is of the Abel form as given in eq. (3.11) of Dejonghe & Laurent (1991). If we substitute this form in (4.52), we obtain, after changing the order of integration and evaluating the integral with respect to $`\lambda `$, again a single Abel integral that is equivalent to the expression for $`T_{\lambda \lambda }`$ that follows from eq. (3.10) of by Dejonghe & Laurent (1991). Using the relations (4.49) and the corresponding subsidiary conditions for $`g_1`$, $`g_2`$ and $`g_3`$, it can be shown that the general triaxial solution (4.45) gives the stress components correctly.

#### 4.6.8 Numerical test

We have numerically implemented the general triaxial solution (4.45), and tested it on a polytrope dynamical model, for which the DF depends only on energy $`E`$ as $`f(E)E^{n3/2}`$, with $`n>\frac{1}{2}`$. Integration of this DF over velocity $`v`$, with $`E=V\frac{1}{2}v^2`$ for a potential $`V0`$, shows that the density $`\rho (V)^n`$ (e.g., Binney & Tremaine 1987, p. 223). This density is not consistent with the Stäckel potentials we use but, as noted in §2.3, the Jeans equations do not require self-consistency. The first velocity moments and the mixed second moments of the DF are all zero. The remaining three moments all equal $`V/(n+1)`$, so that the isotropic stress of the polytrope model $`T_{\mathrm{pol}}(V)^{n+1}`$.

We take the potential $`V`$ to be of Stäckel form $`V_S`$ (2.3), and consider two different choices for $`G(\tau )`$ in (2.4). The first is the simple form $`G(\tau )=GM/(\sqrt{\tau }+\sqrt{\alpha })`$ that is related to Hénon’s isochrone (de Zeeuw & Pfenniger 1988). The second is the form for the perfect ellipsoid, for which $`G(\tau )`$ is given in Z85 in terms of complete elliptic integrals. The partial derivatives of $`V_S(\lambda ,\mu ,\nu )`$, that appear in the weights $`g_1`$, $`g_2`$ and $`g_3`$, can be obtained in terms of $`G(\tau )`$ and its derivative in a straightforward way by using the expressions derived by de Zeeuw et al. (1986).

The calculation of the stresses is done in the following way. We choose the polytrope index $`n`$, and fix the triaxial Stäckel model by choosing $`\alpha `$, $`\beta `$ and $`\gamma `$. This gives $`T_{\mathrm{pol}}`$. Next, we obtain successively the stresses $`T_{\lambda \lambda }`$, $`T_{\mu \mu }`$ and $`T_{\nu \nu }`$ from the general triaxial solution (4.45) by numerical integration, where the relation between $`S_{\tau \tau }`$ and $`T_{\tau \tau }`$ is given by (4.1). We first fix the upper integration limits $`\lambda _e`$, $`\mu _e`$ and $`\nu _e`$. All integrands contain the singular solutions (4.42), that involve the homogeneous solutions $`A`$, $`B`$, $`C`$, $`F`$, $`G`$, $`H`$ and $`I`$, for which we numerically evaluate the single quadratures (eq. \[4.5.1\], \[4.38\] and \[4.5.3\]). The weights $`g_1`$, $`g_2`$ and $`g_3`$ (4.3) involve the polytrope density and Stäckel potential. This leaves the boundary stresses in the integrands, for which we use the polytrope stress $`T_{\mathrm{pol}}`$ that follows from the choice of the DF, evaluated at the corresponding boundary surfaces. We then evaluate the general solution away from these boundaries, and compare it with the known result.

We carried out the numerical calculations for different choices of $`n`$, $`\alpha `$, $`\beta `$ and $`\gamma `$ and at different field points $`(\lambda ,\mu ,\nu )`$. In each case the resulting stresses $`T_{\lambda \lambda }`$, $`T_{\mu \mu }`$ and $`T_{\nu \nu }`$ – independently calculated – were equivalent to high precision and equal to $`T_{\mathrm{pol}}`$. This agreement provides a check on the accuracy of both our formulae and their numerical implementation, and demonstrates the feasibility of using our methods for computing triaxial stress distributions. That will be the subject of a follow-up paper.

## 5 Discussion and conclusions

Eddington (1915) showed that the velocity ellipsoid in a triaxial galaxy with a separable potential of Stäckel form is everywhere aligned with the confocal ellipsoidal coordinate system in which the equations of motion separate. Lynden–Bell (1960) derived the three Jeans equations which relate the three principal stresses to the potential and the density. They constitute a highly-symmetric set of first-order partial differential equations in the three confocal coordinates. Solutions were found for the various two-dimensional limiting cases, but with methods that do not carry over to the general case, which, as a consequence, remained unsolved.

In this paper, we have introduced an alternative solution method, using superposition of singular solutions. We have shown that this approach not only provides an elegant alternative to the standard Riemann–Green method for the two-dimensional limits, but also, unlike the standard methods, can be generalised to solve the three-dimensional system. The resulting solutions contain complete (hyper)elliptic integrals which can be evaluated in a straightforward way. In the derivation, we have recovered (and in some cases corrected) all previously known solutions for the various two-dimensional limiting cases with more symmetry, as well as the two special solutions known for the general case, and have also clarified the restrictions on the boundary values. We have numerically tested our solution on a polytrope model.

The general Jeans solution is not unique, but requires specification of principal stresses at certain boundary surfaces, given a separable triaxial potential, and a triaxial density distribution (not necessarily the one that generates the potential). We have shown that these boundary surfaces can be taken to be the plane containing the long and the short axis of the galaxy, and, more specifically, the part that is crossed by all three families of tube orbits and the box orbits. This is not unexpected, as HZ92 demonstrated that the phase-space distribution functions of these triaxial systems are defined by specifying the population of each of the three tube orbit families in a principal plane. Once the tube orbit populations have been defined in this way, the population of the box orbits is fixed, as it must reproduce the density not contributed by the tubes, and there is only one way to do this. While HZ92 chose to define the population of inner and outer long axis tubes in a part of the $`(x,z)`$-plane, and the short axis tubes in a part of the $`(y,z)`$-plane, it is in fact also possible to specify all three of them in the appropriate parts of the $`(x,z)`$-plane, just as is needed for the stresses.

The set of all Jeans solutions (4.45) contains all the stresses that are associated with the physical distribution functions $`f0`$, but, as in the case of spherical and axisymmetric models, undoubtedly also contains solutions which are unphysical, e.g., those associated with distribution functions that are negative in some parts of phase space. The many examples of the use of spherical and axisymmetric Jeans models in the literature suggest nevertheless that the Jeans solutions can be of significant use.

While triaxial models with a separable potential do not provide an adequate description of the nuclei of galaxies with cusped luminosity profiles and a massive central black hole, they do catch much of the orbital structure at larger radii, and in some cases even provide a good approximation of the galaxy potential. The solutions for the mean streaming motions, i.e., the first velocity moments of the distribution function, are quite helpful in understanding the variety of observed velocity fields in giant elliptical galaxies and constraining their intrinsic shapes (e.g., Statler 1991, 1994b; Arnold et al.1994; Statler, DeJonghe & Smecker-Hane 1999; Statler 2001). We expect that the projected velocity dispersion fields that can be derived from our Jeans solutions will be similarly useful, and, in particular, that they can be used to establish which combinations of viewing directions and intrinsic axis ratios are firmly ruled out by the observations. As some of the projected properties of the Stäckel models can be evaluated by analytic means (Franx 1988), it is possible that this holds even for the intrinsic moments considered here. Work along these lines is in progress.

The solutions presented here constitute a significant step towards completing the analytic description of the properties of the separable triaxial models, whose history by now spans more than a century. It is remarkable that the entire Jeans solution can be written down by means of classical methods. This suggests that similar solutions can be found for the higher dimensional analogues of (2.16), most likely involving hyperelliptic integrals of higher order. It is also likely that the higher-order velocity moments for the separable triaxial models can be found by similar analytic means, but the effort required may become prohibitive.

## acknowledgments

This paper owes much to Donald Lynden–Bell’s enthusiasm and inspiration. This work begun during a sabbatical visit by CH to Leiden Observatory in 1992, supported in part by a Bezoekersbeurs from NWO, and also by NSF through grant DMS 9001404. CH is currently supported by NSF through grant DMS 0104751. This research was supported in part by the Netherlands Research School for Astronomy NOVA. The authors gratefully acknowledge stimulating discussions with Wyn Evans during the initial phases of this work.

## Appendix A Solving for the difference in stress

We compare our solution for the stress components $`T_{\lambda \lambda }`$ and $`T_{\mu \mu }`$ with the result derived by EL89. They combine the two Jeans equations (2.5.3) into the single equation

$$\frac{^2\mathrm{\Delta }}{\lambda \mu }+(\frac{}{\mu }\frac{}{\lambda })\frac{\mathrm{\Delta }}{2(\lambda \mu )}=\frac{\rho }{\lambda }\frac{V_S}{\mu }\frac{\rho }{\mu }\frac{V_S}{\lambda },$$
(A.1)
for the difference $`\mathrm{\Delta }T_{\lambda \lambda }T_{\mu \mu }`$ of the two stress components. Eq. (A.1) is of the form

$$^{}\mathrm{\Delta }=\frac{\rho }{\lambda }\frac{V_S}{\mu }\frac{\rho }{\mu }\frac{V_S}{\lambda },$$
(A.2)
where $`^{}`$ is the adjoint operator defined in eq. (3.6). As in §3.1, eq. (A.1) can be solved via a Riemann–Green function.

### A.1 The Green’s function

In order to obtain the Riemann–Green function $`𝒢^{}`$ for the adjoint operator $`^{}`$, we use the reciprocity relation (Copson 1975, §5.2) to relate it to the Riemann–Green function $`𝒢`$, derived in §3.1.2 for $``$. With $`c_1=c_2=\frac{1}{2}`$ in this case, we get

$`𝒢^{}(\lambda ,\mu ;\lambda _0,\mu _0)`$ $`=`$ $`𝒢(\lambda _0,\mu _0;\lambda ,\mu )`$ (A.3)
$`=`$ $`\left({\displaystyle \frac{\lambda _0\mu _0}{\lambda \mu }}\right)^{\frac{1}{2}}{}_{2}{}^{}F_{1}^{}({\scriptscriptstyle \frac{1}{2}},{\scriptscriptstyle \frac{3}{2}};1;w),`$
where $`w`$ as defined in (3.16). EL89 seek to solve eq. (A.2) using a Green’s function $`G`$ which satisfies the equation

$$^{}G=\delta (\lambda _0\lambda )\delta (\mu _0\mu ).$$
(A.4)
That they impose the same boundary conditions that we do is evident from their remark that, if $`^{}`$ were the simpler operator $`^2/\lambda \mu `$, $`G`$ would be $`(\lambda _0\lambda )(\mu _0\mu )`$. This is the same result as would be obtained by the singular solution method of §3.2, which, as we showed there, is equivalent to the Riemann–Green analysis. Hence their $`G`$ should match the $`𝒢^{}`$ of eq. (A.3). We show in §A.3 that it does not.

### A.2 Laplace transform

We use a Laplace transform to solve (A.4) because the required solution is that to an initial value problem to which Laplace transforms are naturally suited. The PDE is hyperbolic with the lines $`\lambda =\mathrm{const}`$ and $`\mu =\mathrm{const}`$ as characteristics, and its solution is non-zero only in the rectangle bounded by the characteristics $`\lambda =\lambda _0`$ and $`\mu =\mu _0`$, and the physical boundaries $`\lambda =\alpha `$ and $`\mu =\beta `$ (Fig. 7). We introduce new coordinates

$$\xi =(\lambda \mu )/\sqrt{2},\eta =(\lambda +\mu )/\sqrt{2},$$
(A.5)
so that eq. (A.4) simplifies to

$$^{}G\frac{^2G}{\eta ^2}\frac{^2G}{\xi ^2}\frac{}{\xi }\left(\frac{G}{\xi }\right)=2\delta (\xi \xi _0)\delta (\eta \eta _0),$$
(A.6)
where $`\xi _0=(\lambda _0\mu _0)/\sqrt{2}`$ and $`\eta _0=(\lambda _0+\mu _0)/\sqrt{2}`$ are the coordinates of the source point. The factor of 2 arises from the transformation of the derivatives; the product of the delta functions in (A.4) transforms into that of (A.6) because the Jacobian of the transformation (A.5) is unity. The reason for our choice of $`\eta `$ is that $`G0`$ for $`\eta <\eta _0`$, that is $`\lambda +\mu >\lambda _0+\mu _0`$. Hence $`\eta `$ is a time-like variable which increases in the direction in which the non-zero part of the solution propagates. We take a Laplace transform in $`\stackrel{~}{\eta }=\eta \eta _0`$, and transform $`G(\xi ,\eta )`$ to

$$\widehat{G}(\xi ,p)=\underset{0}{\overset{\mathrm{}}{}}e^{p\stackrel{~}{\eta }}G(\xi ,\stackrel{~}{\eta })d\stackrel{~}{\eta }.$$
(A.7)
There are two equally valid ways of taking proper account of the $`\delta (\eta \eta _0)`$ in taking the Laplace transform of equation (A.6). One can either treat it as $`\delta (\stackrel{~}{\eta }0+)`$, in which case it has a Laplace transform of 1, or one can treat it as $`\delta (\stackrel{~}{\eta }0)`$, in which case it contributes a unit initial value to $`G/\eta `$ which must be included in the Laplace transform of $`^2G/\eta ^2`$ (Strauss 1992). Either way leads to a transformed equation for $`\widehat{G}(\xi ,p)`$ of

$$p^2\widehat{G}\frac{\mathrm{d}^2\widehat{G}}{\mathrm{d}\xi ^2}\frac{\mathrm{d}}{\mathrm{d}\xi }\left(\frac{\widehat{G}}{\xi }\right)=2\delta (\xi \xi _0).$$
(A.8)
The homogeneous part of eq. (A.8) is the modified Bessel equation of order one in the variable $`p\xi `$. Two independent solutions are the modified Bessel functions $`I_1`$ and $`K_1`$. The former vanishes at $`\xi =0`$ and the latter decays exponentially as $`\xi \mathrm{}`$. We need $`\widehat{G}`$ to decay exponentially as $`\xi \mathrm{}`$ because $`G(\xi ,\eta )`$ vanishes for $`\stackrel{~}{\eta }<\xi \xi _0`$, and hence its Laplace transform $`\widehat{G}`$ is exponentially small for large $`\xi `$. We also need $`\widehat{G}`$ to vanish at $`\xi =0`$ where $`\lambda =\mu `$. The focus at which $`\lambda =\mu =\alpha `$ is the only physically relevant point at which $`\xi =0`$. It lies on a boundary of the solution region in the $`\lambda _0\alpha `$ limit (Fig. 7). The focus is a point at which the difference $`\mathrm{\Delta }`$ between the stresses vanishes, and hence $`G`$ and $`\widehat{G}`$ should vanish there. The delta function in eq. (A.8) requires that $`\widehat{G}`$ be continuous at $`\xi =\xi _0`$ and that $`\mathrm{d}\widehat{G}/\mathrm{d}\xi `$ decrease discontinuously by 2 as $`\xi `$ increases through $`\xi =\xi _0`$. Combining all these requirements, we obtain the result

$$\widehat{G}(\xi ,p)=\{\begin{array}{cc}2\xi _0K_1(p\xi )I_1(p\xi _0),\hfill & \xi _0\xi <\mathrm{},\hfill \\ 2\xi _0K_1(p\xi _0)I_1(p\xi ),\hfill & 0\xi \xi _0.\hfill \end{array}$$
(A.9)
We use the Wronskian relation $`I_1(x)K_1^{}(x)I_1^{}(x)K_1(x)=1/x`$ (eq. \[9.6.15\] of Abramowitz & Stegun 1965) in calculating the prefactor of the products of modified Bessel functions. The inversion of this transform is obtained from formula (13.39) of Oberhettinger & Badii (1973) which gives

$$G(\xi ,\stackrel{~}{\eta })=\{\begin{array}{cc}\sqrt{\frac{\xi _0}{\xi }}{}_{2}{}^{}F_{1}^{}({\scriptscriptstyle \frac{1}{2}},{\scriptscriptstyle \frac{3}{2}};1;w),\hfill & |\xi _0\xi |\stackrel{~}{\eta }\xi _0+\xi ,\hfill \\ 0,\hfill & \mathrm{}<\stackrel{~}{\eta }<|\xi _0\xi |,\hfill \end{array}$$
(A.10)
we have (cf. eq. \[3.16\])

$$w\frac{\stackrel{~}{\eta }^2(\xi _0\xi )^2}{4\xi _0\xi }=\frac{(\lambda _0\lambda )(\mu _0\mu )}{(\lambda _0\mu _0)(\lambda \mu )}.$$
(A.11)
The second case of eq. (A.10) shows that $`G`$ does indeed vanish outside the shaded sector $`\lambda <\lambda _0`$, $`\mu <\mu _0`$. The first case shows that it agrees with the adjoint Riemann–Green function $`𝒢^{}`$ of (A.3) which was derived from the analysis of §3.1.

### A.3 Comparison with EL89

EL89 use variables $`s=\eta `$ and $`t=\xi `$, whereas we avoided using $`t`$ for the non-time-like variable. They consider the Fourier transform

$$\overline{G}(\xi ,k)=\underset{\mathrm{}}{\overset{\mathrm{}}{}}e^{ik\stackrel{~}{\eta }}G(\xi ,\stackrel{~}{\eta })d\stackrel{~}{\eta }.$$
(A.12)
Because $`G0`$ for $`\stackrel{~}{\eta }0`$, we can rewrite our Laplace transform as their Fourier transform. Setting $`p=ik`$ gives $`\overline{G}(\xi ,k)=i\widehat{G}(\xi ,ik)`$, and using the formulas $`I_1(x)=J_1(ix)`$ and $`K_1(x)=\frac{1}{2}\pi iH_1^{(1)}(ix)`$, eq. (A.9) yields

$$\overline{G}(\xi ,k)=\{\begin{array}{cc}\pi i\xi _0H_1^{(1)}(k\xi )J_1(k\xi _0),\hfill & \xi _0\xi <\mathrm{},\hfill \\ \pi i\xi _0H_1^{(1)}(k\xi _0)J_1(k\xi ),\hfill & 0\xi \xi _0.\hfill \end{array}$$
(A.13)
This formula differs from the solution for the Fourier transform given in eq. (70) of EL89. The major difference is that their solution has Hankel functions of the second kind $`H_1^{(2)}(kt)=H_1^{(2)}(k\xi )`$ where ours has $`J_1`$ Bessel functions. Consequently their solution has an unphysical singularity at $`t=\xi =0`$, and so, in our opinion, is incorrect. Our solution, which was devised to avoid that singularity, gives a result which matches that derived by Riemann’s method in §3.1.

### A.4 The solution for $`\mathrm{\Delta }`$

The solution for $`\mathrm{\Delta }`$ using the adjoint Riemann–Green function is given by eq. (3.14) with $`𝒢`$ replaced by $`𝒢^{}`$ and the sign of $`c_2`$ changed for the adjoint case (Copson 1975). The hypergeometric function of eq. (A.3) for $`𝒢^{}`$ is expressible in terms of complete elliptical integrals as

$${}_{2}{}^{}F_{1}^{}({\scriptscriptstyle \frac{1}{2}},{\scriptscriptstyle \frac{3}{2}};1;w)=\frac{2}{\pi }[E(w)+2wE^{}(w)].$$
(A.14)
Hence, the solution for the difference $`\mathrm{\Delta }`$ between the two principal stresses is given by

$$\begin{array}{c}\mathrm{\Delta }(\lambda ,\mu )=\frac{2}{\pi (\lambda \mu )^{\frac{1}{2}}}\{\hfill \\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}d\lambda _0\underset{\mu }{\overset{\alpha }{}}d\mu _0[E(w)+2wE^{}(w)](\lambda _0\mu _0)^{\frac{1}{2}}(\frac{\rho }{\lambda _0}\frac{V_S}{\mu _0}\frac{\rho }{\mu _0}\frac{V_S}{\lambda _0})\\ \hfill \underset{\lambda }{\overset{\mathrm{}}{}}\mathrm{d}\lambda _0[E(w)+2wE^{}(w)\underset{\mu _0=\alpha }{]}\frac{\mathrm{d}}{\mathrm{d}\lambda _0}[(\lambda _0+\alpha )^{\frac{1}{2}}\mathrm{\Delta }(\lambda _0,\alpha )]\}.\end{array}$$
(A.15)
The determined reader can verify, after some manipulation, that this expression is equivalent to the difference between the separate solutions (3.21a) and (3.21b), derived in §3.1.

### Note added in manuscript

We agree with the amendment to our method of solution for $`\mathrm{\Delta }`$ given in Appendix A.4. Our Green’s function, while solving the differential equation, had the wrong boundary conditions.
N.W. Evans & D. Lynden-Bell