Contents

UTTG-08-12

TCC-010-12

Kinetic Theory of Collisionless Self-Gravitating Gases:

II. Relativistic Corrections in Galactic Dynamics

J. Ramos-Caro111javier@ime.unicamp.br, C. A. Agón222cesar.agon@nucleares.unam.mx and J. F. Pedraza333jpedraza@physics.utexas.edu

Campinas, São Paulo 13083-859, Brazil

Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México,

Apartado Postal 70-543, México D.F. 04510, México

Theory Group, Department of Physics and Texas Cosmology Center,

University of Texas, 1 University Station C1608, Austin, TX 78712, USA

Abstract

In this paper we study the kinetic theory of many-particle astrophysical systems imposing axial symmetry and extending our previous analysis in Phys. Rev. D 83, 123007 (2011). Starting from a Newtonian model describing a collisionless self-gravitating gas, we develop a framework to include systematically the first general relativistic corrections to the matter distribution and gravitational potentials for general stationary systems. Then, we use our method to obtain particular solutions for the case of the Morgan & Morgan disks. The models obtained are fully analytical and correspond to the post-Newtonian generalizations of classical ones. We explore some properties of the models in order to estimate the importance of post-Newtonian corrections and we find that, contrary to the expectations, the main modifications appear far from the galaxy cores. As a by-product of this investigation we derive the corrected version of the tensor virial theorem. For stationary systems we recover the same result as in the Newtonian theory. However, for time dependent backgrounds we find that there is an extra piece that contributes to the variation of the inertia tensor.

## 1 Introduction

The dynamics and evolution of collisionless stellar ensembles is a subject of great interest in astrophysics since they are the primary tool for comparisons of observations and theory in galactic dynamics [1]. Such systems are generally composed by billions of stars, so it is neither practical nor worthwhile to follow the orbit of each particle in the ensemble. Most testable predictions depend only on the distribution function (DF), a quantity that determines the probability of finding a single star in a given phase-space volume around the position and the velocity . The dynamics of the DF follows from the appropriate kinetic equation and it in turn determines the statistical evolution of the system.

In the framework of the general theory of relativity (GR) it is assumed that the DF satisfies the general relativistic version of the Fokker-Planck equation [2, 3, 4] or the collisionless Boltzmann equation (CBE) [5, 6]. The first one is devoted to systems in which local gravitational encounters dominate in their evolution whereas the latter is useful to study systems sufficiently smooth, so that they may be considered to be collisionless [1]. One can actually consider systems in which a number of particle species can collide and produce different species. This is how the formation of the light elements in the big bang nucleosynthesis is calculated (see [7] for a review). However, in systems such as galaxies and galaxy clusters, physical encounters between the stars are very rare, and the effect of gravitational collisions can be neglected for times far longer than the age of the universe. Thus, for these systems, the CBE provides a very good approximation.

For a typical galaxy, the relaxation time, , is arbitrarily large in comparison with the crossing time, . This means that they can be approximated as a continuum rather than concentrated into nearly pointlike stars. Now, it is commonly assumed that the main contribution of the mass in a galaxy is concentrated in an axisymmetric flat distribution [1]. For this reason the obtention of idealized thin disk models has been a problem of great astrophysical relevance. In this case, the most straightforward way to construct a self-consistent model is by means of finding the DF for a system with a known gravitational interactions and matter distribution. Since the mass density of the system is defined by the integration of the DF over the velocity, the problem of finding a DF is that of solving an integral equation (see [8, 9, 10, 11, 12, 13] and the references therein). At present we have at disposal a variety of self-consistent galaxy models: [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25].

Now, even though for most systems under consideration Newtonian gravity is believed to be dominant, general relativistic corrections might play an important role in their evolution. As a matter of fact, in recent years it has been an increasing interest in the incorporation of GR in the description of these systems, and to date we have a variety of fully relativistic galaxy models: [26, 27, 28, 29, 30, 31, 32, 33, 34], among others. Perhaps the principal reason of including GR corrections in galactic dynamics, is the hypothesis that it is possible to overcome the problem of the rotation curves predicted by the Newtonian theory. While, some authors argue that by using GR the inclusion of a dark matter halo is unnecessary at galactic scales (see for instance [35, 36, 37, 38, 39]), several publications have pointed out that this is not entirely true [40, 41, 42, 43, 44]. In particular, the authors of [45] presented a model in which the percentage of dark matter needed to explain flat rotation curves turns out to be less than the required by the Newtonian theory111Strictly speaking, the aim of introducing a dark matter halo is not necessarily to reproduce a flat rotation curve but generate a curve that rises and then declines at particular paces [46].. It is important to point out that currently there are alternative approaches to GR which address the problem of rotation curves in spiral galaxies, as for example the so-called gravity (see [47], for references).

Despite the fact that the relativistic contributions do not solve completely the problem of rotation curves in galaxies, it seems that they do introduce significant corrections. Thus, in order to estimate the effects on the various observables we are interested in, it would be nice to have a framework to include systematically general relativistic corrections to a given Newtonian model. The post-Newtonian approximation is perfectly suited for this purpose. The appropriate scheme that describes the effects of the first corrections beyond the Newtonian theory, was first formulated in [48, 49, 50] (see [51] for a textbook analysis) and it is known as the first post-Newtonian (1PN) approximation. This approach holds if the particles in the system are moving nonrelativistically () (as in the case of a star moving around a typical galaxy) and gives the corrections up to order , where is a typical velocity in the system and is the speed of light. Currently, higher order PN approximations have been developed because of the increasing interest around kinematics and associated emission of gravitational waves by binary pulsars, neutron stars and black holes, with promising candidates for detectors such as LIGO, VIRGO and GEO600 (see [52, 53] for references).

Based in the above considerations, we recently started a general study of self-gravitating gases in the collisionless regime and, as a first step, we derived a version of the CBE that accounts for the first general relativistic corrections [54]. With this tool in hand, we obtained the 1PN version of the Eddington’s polytropes, starting from an ergodic DF proportional to . The purpose of this paper is twofold. First, to implement a similar procedure in the axially symmetric case in order to setup our general framework. And second, to obtain a new set of self-consistent models starting from a Newtonian “seed” and study the impact of relativistic corrections on the various observables.

The rest of the paper is organized as follows. In section 2 we present a brief overview about the basics of the 1PN approximation, revisiting the field equations, as well as the kinetic theory for arbitrary self-gravitating systems. In section 3 we show the fundamental equations defining self-consistent models with post-Newtonian corrections. We start dealing with arbitrary systems but then we focus on discoidal configurations with axial symmetry, in order to prepare the ground to construct 1PN galaxy models in section 4. Finally, we summarize the principal results in section 5.

## 2 General Framework

### 2.1 The 1PN approximation

The post-Newtonian approximation has been reviewed carefully in a number of references (see for example [51]). However, we will include here the basic definitions and relations for completeness.

First off, note that in Newtonian mechanics the typical kinetic energy is roughly of the same order of magnitude as the typical potential energy, and thus

 ¯v2∼ϕ,

where is the mean velocity in the system. The idea is then to express all physical quantities in terms of a series expansion of the small parameter and keep the leading order beyond the Newtonian theory. The first quantity to consider is the spacetime itself: any manifold can be considered to be locally flat so, for particles that are moving nonrelativistically, we proceed to express the metric tensor as

 g00 = −1+\lx@stackrel2g00+\lx@stackrel4g00+⋯, gij = δij+\lx@stackrel2gij+\lx@stackrel4gij+⋯, (1) g0i = \lx@stackrel1g0i+\lx@stackrel3g0i+\lx@stackrel5g0i+⋯,

where the symbol denotes the term in of order . Odd powers of appear in because these components must change sign under time-reversal transformation .

It is natural to assume a similar expansion for the components of the energy momentum tensor. From their interpretation as the energy density, momentum flux, and energy flux, we expect that

 T00 = \lx@stackrel0T00+\lx@stackrel2T00+⋯, Tij = \lx@stackrel2Tij+\lx@stackrel4Tij+⋯, (2) T0i = \lx@stackrel1T0i+\lx@stackrel3T0i+⋯.

These expansions lead to a consistent solution of Einstein field equations.

Working in harmonic coordinates (i.e., coordinates such that ) and to our order of approximation, the various components of the metric tensor can be expressed in terms of the Newtonian potential and post-Newtonian potentials and as

 \lx@stackrel2g00=−2ϕ/c2,\lx@stackrel4g00=−2(ϕ2+ψ)/c4,\lx@stackrel2gij=−2ϕδij/c2,\lx@stackrel1g0i=0,\lx@stackrel3g0i=ξi/c3. (3)

Thus, the Einstein equations reduce to

 ∇2ϕ = 4πG\lx@stackrel0T00, (4) ∇2ψ = 4πGc2(\lx@stackrel2T00+\lx@stackrel2Tii)+∂2ϕ∂t2, (5) ∇2ξi = 16πGc\lx@stackrel1T0i, (6)

along with the coordinate condition

 4∂ϕ∂t+∇⋅\boldmathξ=0. (7)

One can also consider the motion of test particles in a given background. For general potentials , , and one finds that the free falling particle obeys the equation

 dvdt=−∇ϕ−1c2[∇(2ϕ2+ψ)+∂\boldmathξ∂t−v×(∇×\boldmathξ)−3v∂ϕ∂t−4v(v⋅∇ϕ)+v2∇ϕ], (8)

which partially resembles the mathematical structure of the Lorentz force experienced by a charged particle, with velocity , in the presence of an electromagnetic field. Such law of motion will determine, for instance, the rotation curve corresponding to a given galactic model (see Appendix A for details).

It is instructive to point out that the equations of motion (8) can be derived from the Lagrangian [51]

 L=v22−ϕ−1c2(ϕ22+3ϕv22−v48+ψ−v⋅\boldmathξ% ). (9)

For stationary spacetimes, the potentials are independent of time and the associated Hamiltonian , is a conserved quantity that can be interpreted as the 1PN generalization of the classical energy:

 E=v22+ϕ+1c2(3v48−3v2ϕ2+ϕ22+ψ). (10)

Note that this expression is independent of the vector field .

If the source of gravitation is endowed with axial symmetry, the -component of the angular momentum is an additional integral of motion. In cylindrical coordinates we obtain that, for -independent potentials, the quantity

 Lz=Rvφ+1c2[Rvφ(v22−3ϕ)+Rξφ] (11)

can be interpreted as the 1PN generalization of the azimuthal angular momentum. In this case plays a role, through its rotational component.

### 2.2 1PN statistical mechanics

From a statistical point of view, the state of the system can be determined by its DF, , depending on the spatial coordinates, velocity and time. Now, as mentioned in the Introduction, for applications in galactic dynamics it is commonly assumed that the encounters between particles are negligible and hence, the evolution of the stellar system must obey the so-called collisionless Boltzmann equation. In 1PN approximation, such relation can be written as [54]

 ∂F∂t+vi∂F∂xi−∂ϕ∂xi∂F∂vi+1c2(v22−ϕ)(∂F∂t+vi∂F∂xi)+1c2[4vivj∂ϕ∂xj−(3v22+3ϕ)∂ϕ∂xi−vj(∂ξi∂xj−∂ξj∂xi)+3vi∂ϕ∂t−∂ψ∂xi−∂ξi∂t]∂F∂vi=0. (12)

For situations in which encounters play a dominant role, the right-hand side of the above expression can be replaced by a collisional term of the Fokker-Planck type [55].

If the self-gravitating system is in a stationary state, and the DF can be expressed as a function of the energy (a first integral of motion). Moreover, if the system is endowed with additional symmetries we can expect that, as in Newtonian theory, the stationary solutions of (12) are also functions of the remaining integrals of motion. In other words, there is a 1PN version of the Jeans theorem, which simply reflects the fact that equation (12) can be rewritten as [54].

The DF describes completely the sate of the system. We can extract from it all the relevant physical quantities, like the mass density, the mean velocity, the velocity-dispersion tensor, etc. For instance, by taking the moments of equation (12), we could in principle obtain some useful results regarding the hydrodynamics of self-gravitating systems at 1PN order. This was first derived by Chandrasekhar [56] using a different approach so we will refrain from writing out these results here, since they are not particularly illuminating.

Another important fact that can be derived with the help of the DF, is the 1PN virial theorem (see appendix B for details). This theorem shows that there exists a linear relation between the variation of the moment of inertia tensor , the kinetic energy tensor and the potential energy tensor , according to:

 d2Iikdt2=2Kik+Wik+12ddt∫d3x\lx@stackrel0ρc2∂ϕ∂txixk. (13)

This 1PN virial equation has a similar form as the Newtonian one, except for the last term in the right-hand side. The temporal variation of contributes to the variation of . In particular, for stationary states we have

 W=−2K,

as in the Newtonian case. This is, the total gravitational potential energy (a negative quantity), is two times the total kinetic energy.

## 3 The 1PN Self-gravitation Equations

The integration of the DF over the velocity space, leads to the energy-momentum tensor. It in turn determines the gravitational field via the Einstein equations, so it should be possible to derive a relation between and the potentials , , and . Such relation is represented by a set of couple equations which we shall call the 1PN self-gravitation equations. At first we deal with the general case of arbitrary self-gravitating systems and then we shall consider the special case of razor thin disks.

### 3.1 Fundamental equations for arbitrary systems

We begin by considering the following expression for the energy-momentum tensor, valid for any self-gravitating system:

 Tμν(xi,t)=1c∫UμUνU0F(xi,Ui,t)√−gd3U. (14)

Here is the four-velocity, which is related to the three-velocity by (Greek indices run from 0 to 3 and Latin indices from 1 to 3). From the relation we can compute the quantities we need to derive the 1PN approximation of the energy-momentum tensor. After some calculations, we obtain

 U0c=1+v22c2−ϕc2, (15)

and

 d3U=(1+5v22c2−3ϕc2)d3v (16)

Now, taking into account that , we find that

 T00 = ∫γFd3v, Tij = 1c2∫γvivjFd3v, (17) T0i = 1c∫γviFd3v,

where

 γ=1+3v2c2−6ϕc2 (18)

is the measure of the integration over velocities.

According to the above relations, we expect that the DF can be expanded in power series of as

 F=\lx@stackrel0F+\lx@stackrel2F+⋯, (19)

where is the Newtonian contribution to the DF (which is itself a solution of the classical CBE) and is the first post-Newtonian correction. Plugging (19) into (3.1) leads to the different components of the energy-momentum tensor at the orders required by the 1PN approximation:

 \lx@stackrel0T00 = ∫\lx@stackrel0Fd3v, (20) \lx@stackrel2T00 = 3c2∫(v2−2ϕ)\lx@stackrel0Fd3v+∫\lx@stackrel2Fd3v, (21) \lx@stackrel2Tij = 1c2∫vivj\lx@stackrel0Fd3v, (22) \lx@stackrel1T0i = 1c∫vi\lx@stackrel0Fd3v, (23)

along with , as expected. With these results we are ready to write the 1PN self-gravitation equations for general stationary systems:

 ∇2ϕ = 4πG∫\lx@stackrel0Fd3v, (24) ∇2ψ = 8πG∫(2v2−3ϕ)\lx@stackrel0Fd3v (25) +4πGc2∫\lx@stackrel2Fd3v, ∇2ξi = 16πG∫vi\lx@stackrel0Fd3v. (26)

To summarize, we can say that a stellar system characterized by an equilibrium DF satisfying (12) is described by a matter distribution given by (20)-(23), and gravitational interactions determined by the field equations (4)-(6). In order to provide a self-consistent description, the relations (24)-(26) must be satisfied. All of these equations are written as power expansions in the small parameter and, in consequence, we can clearly distinguish between the Newtonian contribution and the post-Newtonian corrections.

### 3.2 The case of stationary razor thin disks

The so-called razor thin disks are of special interest in modeling a number of axisymmetric galaxies. In this case, the DF depends on velocities and positions in the form

 F=f(R,vR,vφ)δ(z)δ(vz), (27)

where is the Dirac delta function and is a reduced phase-space density describing the stellar population placed on the equatorial plane (for a finite thin disk of radius , must vanish for ). Equations (24) and (25) can be written as

 ∇2ϕ = 4πGδ(z)∫\lx@stackrel0fd2v, ∇2ψ = 4πGδ(z)[∫(4v2−6ϕ)\lx@stackrel0fd2v+c2∫\lx@stackrel2fd2v]

where . Equation (26) requires a little more attention. First start with equation (6) (which is equivalent to (26)). The general solution that vanishes at infinity can be written as [51]

 ξi(x)=−4Gc∫\lx@stackrel1T0i(x′)d3x′|x−x′|, (30)

which, in our case, reduces to

 −4G∫∫vi\lx@stackrel0Fd3x′d3v|x−x′|=−4G∫∫vi\lx@stackrel0fδ(z′)d3x′d2v|x−x′|,

where we have computed the integral with respect to (we use the notation , and ). This expression can be massaged into a useful form by taking into account two facts: (i) the relation between the cartesian components , and the cylindrical components , , i.e. and ; (ii) since we are dealing with stationary axisymmetric systems, the DF is an even function of [1] and in consequence . Thus, we can write

 ∫∫vx\lx@stackrel0fδ(z′)d3x′d2v|x−x′|=−sinφ∫∫vφ\lx@stackrel0fδ(z′)d3x′d2v|x−x′|

and

 ∫∫vy\lx@stackrel0fδ(z′)d3x′d2v|x−x′|=cosφ∫∫vφ\lx@stackrel0fδ(z′)d3x′d2v|x−x′|.

Now, by introducing the relations and in (30), we obtain

 ξRcosφ−ξφsinφ=sinφ∫∫4Gvφ\lx@stackrel0fδ(z′)d3x′d2v|x−x′|

and

 ξRsinφ+ξφcosφ=−cosφ∫∫4Gvφ\lx@stackrel0fδ(z′)d3x′d2v|x−x′|.

Since we assume that is -independent, each of these expressions leads us to the conclusion that

 ξR=0,for finite distributions, (31)

and that is solution of the following equation:

 ∇2ξφ=16πGδ(z)∫vφ\lx@stackrel0fd2v. (32)

The equation for the component can be obtained easily by replacing (27) in (26), and the result is the Laplace equation, Its solution can be determined through condition (7). A straightforward calculation leads to

 ξz(R)=ξzoln(R/Ro), (33)

where and are constants of integration. In the case of distributions with finite extent, we demand as a boundary condition that . In consequence, we have to choose , and hence

 ξz=0,for finite thin disks. (34)

On the other hand, we expect that obeys the collisionless Boltzmann equation in the three-dimensional phase-space . In fact, by introducing (27) in (12) and performing an integration on and , it follows that the distribution obeys the relation

 −(1+v22c2−ϕc2−4Rc2∂ϕ∂R−Rc2vφ∂ξφ∂R)vRvφR∂f∂vφ+(1+v22c2−ϕc2)vR∂f∂R (35) +[(1+v22c2−ϕc2)v2φR(1+3v2φ−5v2R2c2+3ϕc2)∂ϕ∂R−1c2∂ψ∂R+vφc2∂ξφ∂R]∂f∂vR = 0,

where and , and are evaluated at . The above relation is the 1PN version of the Boltzmann equation for an axisymmetric two-dimensional shell, located at the equatorial plane, and in a stationary state. Of course, plays the role of the reduced DF describing the diskoidal shell.

It is straightforward to show that and , given by equations (10) and (11), are solutions of (35). This means that, for axially symmetric systems, any depending on and is solution of the CBE; conversely, any solution of the CBE can always be expressed as a function of and . Thus, any two-integral DF, , provides a complete statistical description for the (two-degree-of-freedom) stellar system. This fact will be very useful for the formulation of post-Newtonian models in the next section.

## 4 Analytical Models for Axisymmetric Galaxies

The purpose of this section is to show how to implement the formalism developed above in order to obtain axially symmetric galaxy models. For the applications we want to consider here we have to take into account further considerations. First of all, recall that in [54] we proved that Jeans theorem remains valid at 1PN order. This means that any equilibrium solution of the CBE depends only on the integrals of motion of the system, and that any function of the integrals yields an equilibrium solution of the CBE. Thus, for stationary systems with axial symmetry, we can restrict ourselves to DFs depending on the energy (10) and the angular momentum (11), which are themselves integrals of (12).

The next step would be to implement the previous restrictions starting from a given Newtonian potential-density pair with a known DF, as was done in [54] for the spherically symmetric case. As a result, one expects two coupled self-gravitation equations, providing a method to determine, from a Newtonian model, its associated post-Newtonian corrections. In practice, the present formalism leads to a two coupled ordinary differential equations in the spherically symmetric case. In the axially symmetric situation however, one ends up with two coupled elliptic partial differential equations (for general volumetric matter distributions).

In this case such equations are much more involved than the ones corresponding to the spherically symmetric case, but the configurations we shall deal here permit us to introduce some additional assumptions to simplify the problem. In the next section, we shall show that a dramatic simplification can be achieved by the consideration of thin discoidal distributions in spheroidal oblate coordinates: instead of getting differential equations, in this case the post-Newtonian corrections can be obtained from simple algebraic equations.

We then present a particular application where the resulting equations can be solved analytically, which means that it is possible to obtain 1PN exact solutions. The importance of these solutions will be evaluated by a comparison between density profiles and rotation curves described by Newtonian theory and the ones predicted by the 1PN approximation. Although focus on the particular models introduced in [19] (revisited by [24]), our framework can be applied to a wider variety of models. In general, the method can be used for situations in which the potentials are separable functions of the spheroidal oblate coordinates.

### 4.1 Hunter’s method in the 1PN approximation

One can find in the literature a number of self-consistent stellar models representing razor thin disks; here we will deal with models belonging to the family of Morgan & Morgan disks [19, 24]. In the Newtonian formulation, they can be obtained by a formalism developed by Hunter [57]. Such procedure (known as the Hunter’s method) provides the surface density of the disks, the gravitational potential, and the circular velocity as a series of elementary functions, by superposing solutions of the Laplace equation in oblate spheroidal coordinates. Hunter’s method can also be implemented in the context of the 1PN approximation as follows.

To begin with, note that in vacuum, the field equations (4)-(6) reduce to three Laplace equations for , , and (remember that for distributions of finite extent). Without loss of generality, we assume that the disk is on the equatorial plane, so we have to impose that the gravitational potentials have symmetry of reflection with respect to the plane , i.e., , , and . Then, it follows that

 ∂ϕ∂z(R,−z)=−∂ϕ∂z(R,z), (36) ∂ψ∂z(R,−z)=−∂ψ∂z(R,z), (37) ∂ξφ∂z(R,−z)=−∂ξφ∂z(R,z), (38)

in agreement with the attractive character of gravitation. We also assume that , , and do not vanish in the disk’s zone, in order to have the corresponding thin distribution of energy-momentum. Such distribution, restricted to a region in the plane (from here on, will denote the disk radius), will be described by a “shell-like” energy-momentum tensor. If we define

 \lx@stackrel0T00 = Σ(R)δ(z), (39) \lx@stackrel2T00+\lx@stackrel2Tii = 1c2σ(R)δ(z), (40) \lx@stackrel1T0φ = 1cΔ(R)δ(z), (41)

for , it follows from Gauss’s Law that

 Σ(R)=12πG(∂ϕ∂z)z=0+, (42) σ(R)=12πG(∂ψ∂z)z=0+, (43) Δ(R)=18πG(∂ξφ∂z)z=0+. (44)

Note that represents the surface mass density of the Newtonian theory (i.e. without relativistic corrections), while plays the role of the surface density of -momentum. On the other hand, is associated both to the pressure and the relativistic corrections to the mass surface density.

The above relations mean that, in order to have a distribution of matter as the described by (42)-(44), we have to demand that

 ∂ϕ∂z(R,0+)≠0,R≤a, (45) ∂ϕ∂z(R,0+)=0,R>a, (46)

with the same requirement for and . At this point it is convenient to introduce oblate spheroidal coordinates, a system that adapts in a natural way to the geometry of the problem. They are related to the cylindrical ones through

 R=a√(1+ζ2)(1−η2), (47) z=aζη, (48)

where and . Note that (i) the disk itself has coordinates , ; (ii) conditions (45)-(46) become

 ∂ϕ∂ζ∣∣∣ζ=0=H(η), (49) ∂ϕ∂η∣∣∣η=0=0, (50)

where is an even function of . The general solution of Laplace’s equation satisfying the above conditions can be written as

 ϕ(ζ,η)=−∞∑n=0A2nq2n(ζ)P2n(η), (51)

where are arbitrary constants, and are the usual Legendre polynomials and the Legendre functions of second kind, respectively. The post-Newtonian potentials and have the same form,

 ψ(ζ,η)=−∞∑n=0B2nq2n(ζ)P2n(η), (52)
 ξφ(ζ,η)=−∞∑n=0C2nq2n(ζ)P2n(η), (53)

but here we have denoted the expansion constants as and . We can derive explicit formulas for , , and in oblate spheroidal coordinates, by introducing (51)-(53) in (42)-(44):

 Σ=12πaGη∗∞∑n=0A2n(2n+1)q2n+1(0)P2n(η∗), (54) σ=12πaGη∗∞∑n=0B2n(2n+1)q2n+1(0)P2n(η∗), (55) Δ=18πaGη∗∞∑n=0C2n(2n+1)q2n+1(0)P2n(η∗), (56)

where represents the value of coordinate inside the disk:

 η∗=√1−R2a2 (57)

Now that we have stated the fundamental structure of models with 1PN corrections, the next step is to demand that the models obtained are self-consistent, i.e., that they have an analytical equilibrium DF that is related consistently to the surface mass distribution. In other words, we have to formulate the corresponding 1PN self-gravitating equations:

 ∞∑n=0~A2nP2n(η∗)η∗ = ∫\lx@stackrel0fd2v, ∞∑n=0~B2nP2n(η∗)η∗ = ∫(4v2−6ϕ)\lx@stackrel0fd2v+c2∫\lx@stackrel2fd2v, (59) ∞∑n=0~C2nP2n(η∗)η∗ = ∫vφ\lx@stackrel0fd2v, (60)

where, for the sake of simplicity, we have defined

 ~A2n=(2n+1)q2n+1(0)2πaGA2n, (61)

the same for , and

 ~C2n=(2n+1)q2n+1(0)8πaGC2n. (62)

Note that, part of the right-hand side of equations (4.1)-(60) is characterized by three fundamental quantities in Newtonian dynamics of self-gravitating systems: (i) the mass surface density, ; (ii) the mean square velocity, , and (iii) the mean circular velocity, . In particular, equation (4.1) is the Newtonian self-gravitation equation.

There is a variety of cases that can be addressed by the formalism presented above. They correspond to situations in which the gravitational fields (i) have symmetry of reflection with respect to the equatorial plane, (ii) are separable functions of the spheroidal oblate coordinates and (iii) correspond to axially symmetric distributions with finite extent. Note that any function with these features can be expressed in the form (51).

### 4.2 The Morgan & Morgan disks with 1PN corrections

In Newtonian gravity, the Generalized Kalnajs Disks [24] are finite thin distributions of matter with surface mass density

 Σm(R)=(2m+1)M2πa2(1−R2a2)m−1/2,m∈N (63)

and gravitational potential , given by (51), with the following expansion constants

 A(m)2n=MGa√π2−2m−1(4n+1)(2m+1)!(2n+1)(m−n)!Γ(m+n+32)q2n+1(0). (64)

These solutions were first obtained by Morgan & Morgan (MM) [19], by solving a boundary value problem in the context of GR. Even so, we refer to the above solutions as MM disks for historical reasons, but keeping in mind that we are dealing with Newtonian gravity.

All the members of this family have interesting features: a monotonically decreasing mass density and a Keplerian rotation curve. The case , corresponding to the well-known Kalnajs disk [20] (see also [1]), is an exception because it describes a self-gravitating disk which rotates as a rigid body.

The superposition of members belonging to the MM family, has been used to obtain new models with more realistic properties. For example, in reference [25], the authors constructed a family of models with approximately flat rotation curves considering a particular combination of MM disks. Another example is the family of models introduced by Letelier in [58], representing flat rings (see also [59] for astrophysical applications). Additionally, in reference [60] it was shown that it is possible to construct models with realistic rotational curves obeying simple polynomial expressions. In particular, the authors constructed models for a number of galaxies of the Ursa Major cluster, and as an application they estimated their corresponding mass distributions. This fact suggest that there exist superpositions of MM members leading to galactic models in agreement with the so-called maximum disk hypothesis [1]. It would be interesting to have in hand the 1PN version of the MM family, along with all of its features.

Here we illustrate how to obtain the post-Newtonian version of the MM disks. We focus in the model, which is characterized by a density

 Σ=5M2πa2η3∗, (65)

but in principle, the same procedure can be applied to any of the remaining members (or linear combinations). From here on, we will drop the subindex of , for simplicity.

It can be shown that the above surface density can be obtained from the DF [61]

 f(E,Lz)=k(ΩLz−E−5Ω2a2/4)−1/4, (66)

where

 k=2√3[10Mπ11a5G3]1/4,Ω=√15πGM32a3. (67)

Note that the DF defined by (66) is function of the Jacobi’s integral, , which can be interpreted as the energy measured from a frame that rotates with constant angular speed (see Appendix C).

In order to obtain a 1PN version of the model, we start from the DF given by (66) but using the post-Newtonian expressions for and , i.e., equations (10) and (11). Thus, at 1PN order we can write

 f=\lx@stackrel0f+\lx@stackrel2f, (68)

with

 \lx@stackrel0f=kJ−1/40and\lx@stackrel2f=−k4J2J−5/40, (69)

where

 J0 = −ϕ−v2/2+ΩRvφ−5Ω2a2/4, J2 = −c−2[ϕ22+ψ−3v2ϕ2+3v48−ΩRvφ(v22−3ϕ)−ΩRξφ]. (71)

These relations determine the 1PN self-gravitation equations through (4.1)-(60) and the integrals

 ∫\lx@stackrel0fd2v = 5M2πa2η3∗, ∫v2\lx@stackrel0fd2v = 75GM2448a3(7−7η2∗+6η4∗)η3∗ ∫vφ\lx@stackrel0fd2v = 58√15GM32πa5√(1−η2∗)η3∗ ∫\lx@stackrel2fd2v = 1c2η∗[GM2a34∑j=0b2jη2j∗+√160M3π3a3G√(1−η2∗)ξφ+32ψ32aGπ2] (75)

where

 b0=−23251024,b2=75512,b4