Phase space mixing in the equatorial plane of a Kerr black hole
Paola Rioseco Olivier Sarbach







^1Instituto de Física y Matemáticas,
Universidad Michoacana de San Nicolás de Hidalgo,
Edificio C-3, Ciudad Universitaria, 58040 Morelia, Michoacán, México.


It is shown that a collisionless, relativistic kinetic gas configuration propagating in the equatorial plane of a Kerr black hole undergoes a relaxation process and eventually settles down to a stationary, axisymmetric configuration surrounding the black hole. The underlying mechanism for this relaxation process is due to phase space mixing, which implies that although the one-particle distribution function f satisfying the collisionless Boltzmann equation is quasi-periodic in time, the associated macroscopic observables computed from averages over f possess well-defined limits as time goes to infinity. The final state of the gas is described by an effective distribution function depending only upon the constants of motion, and it can be determined by an appropriate average of the initial distribution function.





04.20.-q,04.40.-g, 05.20.Dd










Phase space mixing in the equatorial plane of a Kerr black hole
    Paola Rioseco and Olivier Sarbach
    December 30, 2023
===============================================================







§ INTRODUCTION



Phase space mixing plays an important role in a wide range of areas in physics, including galactic dynamics, plasma physics and quantum physics. Roughly speaking, this phenomenon can be understood as the relaxation of the observables associated with a distribution function which is transported along an anharmonic Hamiltonian flow and hence spreads all over phase space due to different orbits having different angular frequencies.

In the context of galactic dynamics, groundbreaking work by Lynden-Bell <cit.> has shown that phase space mixing might be one of the fundamental mechanisms responsible for driving the one-particle distribution function describing the stellar distribution in a galaxy to an equilibrium configuration, although collisions or the exchange of energy between individual stars are negligible. For further important discussions on these topics, see for instance Refs. <cit.>.

In plasma physics, phase mixing has been argued to be the driving mechanism for Landau damping, the relaxation of a charged collisionless gas to a homogeneous configuration. Recently, this explanation has been put on a rigorous basis by the work of Mouhot and Villani <cit.> who established that for the case of a finite box with periodic boundary conditions the known phase-mixing property for the linearized Vlasov-Poisson also occurs in the full Vlasov-Poisson system, without linearization. For generalization of this work to the special-relativistic setting, see <cit.>.

At the quantum level, phase space mixing has recently been applied <cit.> to quenches in Bose-Einstein condensates in order to understand the behavior of the condensate in the vicinity of the saddle point in a double-well potential. Again, it was found that the system relaxes to a steady state due to phase-space mixing. For field-theoretical applications of the mixing phenomenon, see for example <cit.>.

In the present work, we analyze the effects of phase mixing and its associated relaxation process in a general relativistic scenario. More specifically, we consider a collisionless, relativistic gas configuration that is trapped by the gravitational field of a rotating black hole. We restrict ourselves to the simplest case in which the gas configuration is sufficiently thin such that its self-gravity can be neglected and in which the gas is confined to the equatorial plane of the black hole, leaving the discussion of more realistic configurations to future work <cit.>. As a consequence of our assumptions, each individual gas particle follows a bound geodesic trajectory in the equatorial plane of a Kerr black hole background, and an explicit solution representation for the one-particle distribution function can be obtained by representing the geodesic flow in terms of action-angle-like variables, see Section <ref>.

Based on this solution representation, we compute the particle current density four-vector and provide some examples in Section <ref> showing that the particle density measured by a stationary observer outside the event horizon, although fluctuating in time, undergoes damped oscillations and eventually settles down to a constant value. An intuitive explanation for this convergence is given by exhibiting snapshots of the distribution function in the momentum space of the observer at different times which clearly illustrate the mixing in phase space.

To provide a precise mathematical formulation for the mixing property, in Section <ref> we consider an observable N[φ] which is obtained by integrating the one-particle distribution function over a given test function φ on relativistic phase space. The time evolution of this observable is obtained by transporting the test function along the vector field generating the time translation symmetry. Denoting this transported test function by φ_t, we formulate and prove a theorem which shows that, provided a certain determinant condition holds on the support of φ, N[φ_t] converges for t→∞ to the same observable N[φ] with the distribution function replaced by its average over the angle variables. Therefore, apart from providing a rigorous formulation for the mixing property, our theorem allows to predict the final state of the gas configuration by considering the average of the initial distribution function. Finally, we show that the determinant condition is satisfied almost everywhere in phase space, and discuss the applicability of our theorem to the examples in Section <ref>.

In Section <ref> we present the conclusions and main implications of our results, and also give an outlook to future work. Technical details and relevant analytic expressions required for this article are listed in an appendix. Throughout this work we use geometrized units in which the gravitational constant and the speed of light are one.






§ COLLISIONLESS DISTRIBUTION FUNCTION IN TERMS OF ACTION-ANGLE VARIABLES




The geodesic flow describing the motion of free falling, massive particles following (spatially) bound trajectories on a Kerr spacetime may be represented analytically in terms of action-angle-like variables <cit.>. For the purpose of this work, it is sufficient to consider the 3-dimensional spacetime ( M,g) describing the induced geometry on the equatorial plane of a Kerr black hole exterior of mass M > 0 and rotational parameter a satisfying |a|≤ M. In terms of Boyer-Lindquist coordinates (t,r,φ), the metric has the following representation:

    g = -dt^2 + 2M/r( dt - a dφ)^2 + (r^2 + a^2) dφ^2 + r^2/Δ dr^2,
         r > r_H

with Δ := r^2 - 2Mr + a^2  and r_H := M + √(M^2 - a^2) the radius of the event horizon. Due to the time-translation and axial symmetry of ( M,g), the geodesic equations form an integrable Hamiltonian system which is characterized by the free-particle Hamiltonian H and the constants of motion E and L, given by the following functions on the co-tangent  bundle T^* M associated with M:

    H(x,p) := 1/2 g^μν(x) p_μ p_ν,   E(x,p) := -p_t,   L(x,p) := p_φ,     (x,p)∈ T^* M,

which Poisson-commute among themselves. The orbits are confined to the invariant subsets 

    Γ_m,E,L := { (x,p)∈ T^* M :  H(x,p) = -m^2/2,  E(x,p) = E,  L(x,p) = L },

with m > 0. Since only bound orbits are considered, the angular momentum L has to be large enough in magnitude such that L^2 > L_ms^2, with L_ms the angular momentum corresponding to the marginally stable circular orbit <cit.>, and the energy E has to lie inside a certain interval E_min(L) < E < E_max(L) with E_min(L) the energy of the stable circular orbit with angular momentum L and E_max(L)≤ m the maximum of the potential well. In this range, it can be verified that the invariant sets Γ_m,E,L are smooth 3-dimensional submanifolds of T^* M having topology × S^1× S^1. For the following, we focus on the phase space of bounded trajectories Γ_bound, the union of all these invariant submanifolds. 

Using standard tools from Hamiltonian mechanics <cit.> one can introduce action-angle variables (Q^α,J_α) on Γ_bound, see for instance <cit.>. The action variables are defined as[Our definition of J_0 differs from the definition of J_t in <cit.> by a minus sign.]

    J_0 := -1/T∫_0^T p_t dt =  E,    
    J_1 := 1/2π∮ p_φ dφ =  L,    
    J_2 := 1/2π∮ p_r dr = 1/π∫_r_1^r_2√(R(r))/Δ dr,

where here the first integral defining J_0 is performed along an integral curve of the vector field ∂_t (t being the parameter along this curve), the second integral defining J_1 is similarly performed over the closed curve along the vector field ∂_φ,[More precisely, these integral curves are defined with respect to the complete lifts of the Killing vector fields ∂_t and ∂_φ on T^* M, see for example Refs. <cit.> for the most important properties of the complete lift in the context of relativistic kinetic theory.] and the third integral over the closed curve in the (r,p_r)-plane described by the radial equation (Δ p_r)^2 = R(r) with

    R(r) := (E r^2 - aL̂)^2 - Δ(m^2 r^2 + L̂^2),    L̂ := L - a E,

and with r_1 < r_2 the turning points. The angle variables Q^α := ∂ S/∂ J_α are obtained from the generating function

    S(x; J_0,J_1,J_2) = - E t + Lφ + ∫_(r,p_r) p_r dr,

where the integral on the right-hand side should be interpreted as a line integral along the curve (Δ p_r)^2 = R(r) connecting the reference point (r_1,0) to the given point (r,p_r) on this curve. S, Q^α, and J_α admit explicit representations in terms of standard elliptic integrals (details of the derivation will be presented elsewhere <cit.>). In order to describe the result we denote the roots of the fourth-order polynomial R(r) by 0 < r_0 < r_1 < r_2 and introduce the dimensionless variables

    α := a/M,    ε := E/m,    λ := L/M m,    ξ_j := r_j/M,     j=0,1,2.




The curve (Δ p_r)^2 = R(r) is parametrized by the π-periodic angle coordinate χ (see Fig. <ref>) defined by

    r/M   =   ξ_0 + ξ_1 - ξ_0/1 - b^2sin^2χ,
    Δ p_r/M^2 m   =   1/2√(1 - ε^2)(ξ_2 - ξ_1)(ξ_1 - ξ_0)√(ξ_1/ξ_2 - ξ_0)√(1 - k^2sin^2χ)/(1 - b^2sin^2χ)^2sin(2χ),

with b := √((ξ_2 - ξ_1)/(ξ_2 - ξ_0)) and k := √(ξ_0/ξ_1) b, such that 0 < k < b < 1. In terms of these quantities one finds

    J_0 = mε,     
     J_1 = M mλ,     
     J_2 = M m/π[ (1-ε^2)_0 + ε_3 ],
    
           Q^0 = -t + M_0_2(χ) - _2_0(χ)/_0,    
    Q^1 = φ + _0_1(χ) - _1_0(χ)/_0,    
    Q^2 = π_0(χ)/_0,

where the functions _j(χ) and corresponding constants _j := _j(π/2) are defined in Appendix <ref>. The coordinates (Q^α,J_α) provide new symplectic coordinates on Γ_bound, the submanifold of T^* M corresponding to bound orbits. The action variables J_α label the invariant submanifolds Γ_m,E,L, the coordinates (Q^1,Q^2) are 2π-periodic functions providing angles on each S^1-factor of Γ_m,E,L = × S^1× S^1, while Q^0 parametrizes its -factor.

In terms of these action-angle-like variables, the collisionless Boltzmann equation { H, f } = 0 for the one-particle distribution function f assumes the simple form

    ( Ω^0∂/∂ Q^0 + Ω^1∂/∂ Q^1
     + Ω^2∂/∂ Q^2) f = 0,    Ω^α := ∂ H/∂ J_α.

Since H and Ω^α only depend on the action variables J_α, the most general solution of Eq. (<ref>) has the form

    f(x,p) = F(Q^1 - ω^1 Q^0, Q^2 - ω^2 Q^0, J_0,J_1,J_2),

with F(Q^1,Q^2,J_0,J_1,J_2) an arbitrary function which is 2π-periodic in the Q variables and decays sufficiently fast in the J variables, such that the integrals defining the spacetime observables are well-defined. Here the frequencies ω^1 and ω^2 are defined by

    ω^1 := Ω^1/Ω^0 
     = 1/M_1/_2 - ε_0,    ω^2 := Ω^2/Ω^0 = -π/_1ω^1.

Eqs. (<ref>,<ref>,<ref>) provide an explicit solution representation for the distribution function f which will enable us to study the dynamical behavior of generic, time-dependent gas configurations in the next two sections. Note that according to Eq. (<ref>) the distribution function f is axisymmetric if and only if F is independent of Q^1, and it is stationary and axisymmetric if and only if F is independent of Q^1 and Q^2 (in which case it depends only on the action variables J_α). In the remainder of this work, we demonstrate that any distribution function relaxes in time (in some sense made precise in Section <ref>) to such a stationary and axisymmetric configuration.





§ RELAXATION OF SPACETIME OBSERVABLES AND PHASE SPACE MIXING




In this section we analyze the behavior of spacetime observables along the world lines of observers located outside the event horizon. For simplicity and the sake of illustration, in this section we assume that the black hole is non-rotating (although we will come back to the rotating case in the next section). Furthermore, we focus on the particle current density four-vector (see for instance <cit.>)

    J^μ(x) = ∫ f(x,p) p^μ_x(p),    _x(p) := √(-(g^μν(x))) dp_t dp_r dp_ϑ dp_φ,

measured by a static observer in the equatorial plane with fixed spatial coordinates (r,ϑ,φ) = (r_obs,π/2,φ_obs), r_obs > 2M, with respect to its rest frame

    e_0 = 1/√(1-2M/r)∂/∂ t,    
    e_1 = √(1-2M/r)∂/∂ r,    
    e_2 = 1/r∂/∂ϑ,    
    e_3 = 1/r∂/∂φ.

In order to perform the integral over the momentum we first re-express the volume form _x(p) in terms of the conserved quantities m,E,ℓ and ℓ_z, where ℓ and ℓ_z are, respectively, the total and the azimuthal component of the angular momentum. For points x located in the equatorial plane, one obtains

    _x(p) = 1/r^2 dE dp_r dp_ϑ d ℓ _z 
     = dE (m dm) (ℓ dℓ) dσ/√(R(r)),

where we have defined σ by sinσ = ℓ_z/ℓ and the function R(r) is given in Eq. (<ref>). Assuming a kinetic gas distribution of identical particles of positive rest mass m which are confined to the equatorial plane (such that σ = ±π/2 and ℓ_z = L), one obtains the following orthonormal components of the current density:

    J^α := . m^3∫[ f(x,p_+) p_+^α + f(x,p_-) p_-^α] |λ| dε dλ/√((1-ε^2)ξ(ξ - ξ_0)(ξ - ξ_1)(ξ_2 - ξ))|_ξ = r_obs/M,

with p^0_± := mε/√(1 - 2/ξ), p^1_± := √(1 - 2/ξ) p_r±, p^2_± := 0 and p^3_± := mλ/ξ, and where p_± denote the two possible values for the four-momentum p = p^α e_α corresponding to the two solutions of equation (Δ p_r)^2 = R(r). To make further progress, we use the relations (<ref>) to re-express the integral in terms of the turning points ξ_1 and ξ_2. Instead of ξ_1 and ξ_2, it is convenient to use the generalized Keplerian variables (P,e) as in Ref. <cit.>, which are defined by

    ξ_1 = P/1 + e,    ξ_2 = P/1 - e,

such that ξ_0 = 2P/(P-4). Here, the eccentricity e is restricted to the interval 0 < e < 1, the limit e→ 0 representing circular trajectories, and the semi-latus rectum P is restricted to P > 6 + 2e, the limits P→∞ and P→ 6 + 2e (with fixed e) corresponding, respectively, to the Newtonian limit and the innermost stable orbits (ISO) which separate the bound orbits from those that plunge into the black hole. Note that there is a one-to-two correspondence between the parameters (P,e) and the constants of motion (ε,±λ). Note also that given ξ_obs = r_obs/M, the parameters (P,e) are restricted by the conditions ξ_0 < ξ_1 < ξ_obs < ξ_2, which yield

    P_min(e) := max{ 6 + 2e, (1-e)ξ_obs} < P < (1+e)ξ_obs =: P_max(e).

Based on these observations, Eq. (<ref>) can be rewritten as the sum J^α =  J_λ > 0^α +  J_λ <  0^α over the contributions corresponding to gas particles with positive/negative angular momentum λ, where

    J_λ > 0^α = . m^3/2∑_±∫_0^1 de ∫_P_min(e)^P_max(e) dP
     F(Q_±^1 - ω^1 Q_±^0,Q_±^2 - ω^2 Q_±^0,J_0,J_1,J_2) p^α_±/√((1-ε^2)ξ(ξ - ξ_0)(ξ - ξ_1)(ξ_2 - ξ))e√(P)[ (P-6)^2 - 4e^2]/√((P-3-e^2)^5[(P-2)^2 - 4e^2])|_ξ = ξ_obs.

Here, it is understood that all the relevant quantities are expressed in terms of (P,e) which determine the locations of the roots ξ_0,ξ_1,ξ_2 and the conserved quantities (ε,λ) taking the positive sign of λ (and similarly for J_λ < 0^0 where one takes the negative sign of λ). The sum over the ± signs refers to the two possible choices for p_± which, in turn correspond to the two possible values for the angle variables Q_±^α. Note that Q_+^α + Q_-^α = 2π.

In Figs. <ref>, <ref> and <ref> we display the time behavior of the observer's particle density n :=  J^0_λ > 0 for the case of an initial distribution function of the form F(q^1,q^2,J_0,J_1,J_2) = F_Q(q^1,q^2) F_J(P,e), where

    F_Q(q^1,q^2) := 1/Nexp[ -cos^2(q^1)/η_1^2 - cos^2(q^2)/η_2^2],    
    F_J(P,e) := exp[ -(P - 2e - P_0 + 2e_0)^2/(Δ P)^2 - (e - e_0)^2/(Δ e)^2],

N is a positive normalization factor chosen such that ∫_0^2π dq^1 ∫_0^2π dq^2 F_Q(q^1,q^2) = 1, and η_1,η_2,P_0,e_0,Δ P,Δ e are positive constants whose values are given in Table <ref>. Notice that the function F_Q defined in Eq. (<ref>) satisfies F_Q(q^1,q^2) = F_Q(2π - q^1,2π - q^2), which implies that both terms in the sum ∑_± yield the same contribution for J^α_λ > 0 when α=0,3, while J^1_λ > 0 = 0. The observer is located at either (r_obs,φ_obs) = (6M,0) or (9M,0).



Even though the distribution function oscillates in time and does not have a (pointwise) limit, the results in the plots indicate that the observable n converges to a finite value (denoted by n_∞) as t→∞. The theorem in the next section provides the means for determining this asymptotic value. For the moment we note that in case A the convergence is rather fast (exponential decay during the initial period, with relative fluctuations below 0.001 after times larger than 300M), while in the remaining cases the convergence appears to be slower (inverse power-law decay with relative fluctuations below 0.001 after a few thousand M). We have also computed the remaining nontrival component J^3_λ > 0 and the contributions J^α_λ < 0 belonging to negative angular momentum (assuming the same distribution F_J(P,e)) and verified that they exhibit a similar time behavior, with J^3 :=  J^3_λ > 0 +  J^3_λ < 0→ 0.










As shown in the next section, the relaxation process displayed in Figs. <ref>-<ref> is due to phase space mixing. To get an intuitive idea about this phenomenon, in Fig. <ref> we show snapshots of the function F_Q(Q^1 -ω^1 Q^0,Q^2 - ω^2 Q^0) for an observer located at r_obs = 6M and φ_obs = 0 at different times. As is visible from these plots, the geodesic flow (which is volume preserving according to Liouville's theorem) stretches the phase space elements and spreads them over large regions in phase space, inducing the mixing property. As a consequence, averaged (macroscopic) quantities computed from the distribution function, such as the components of the current density in Eq. (<ref>), have the form of an integral over a smooth function multiplied by an oscillating function whose frequency increases unboundedly in time. In the limit t→∞ these oscillations cancel out, and hence one can replace the distribution function with its non-oscillatory part, that is, its average over the angle variables:

    F(Q^1 - ω^1 Q^0, Q^2 - ω^2 Q^0, J_0,J_1,J_2) ⇀F(J_0,J_1,J_2) := 1/(2π)^2∫_0^2π∫_0^2π F(Q^1,Q^2, J_0,J_1,J_2) dQ^1 dQ^2

The precise sense in which this convergence is valid will be explained in the next section. 








 








§ MATHEMATICAL FORMULATION OF THE MIXING PROPERTY




After discussing the intuitive picture behind the mixing phenomenon and the corresponding relaxation process of the spacetime observables, in this section we provide a precise mathematical formulation of this effect which provides a rigorous explanation for the convergence of certain macroscopic observables. Here, a macroscopic observable is described by a test function φ∈ C_0^∞(Γ_bound) on the relativistic phase space Γ_bound, and its associated value is defined by the quantity

    N[φ] := ∫_Γ_bound f(x,p) φ(x,p) _Γ,

with _Γ = dt dφ dr dp_t dp_φ dp_r the canonical volume form on Γ_bound and f(x,p) the one-particle distribution function describing the state of the kinetic gas. Note that our definition (<ref>) is based on a fully covariant (i.e. independent of any choice of foliation or local coordinates) spacetime point of view which includes a time integral, such that N[φ] does not  dependent on time. Therefore, in order to understand the dynamical behavior, one needs to perform a translation of the test function φ along a time direction. Due to the many-fingered nature of time in general relativity there are no such preferred time directions in a general situation. In our case, the presence of the Killing vector field ∂_t induces a natural vector field k̂ on phase space Γ_bound defined as the complete lift of ∂_t (see Refs. <cit.> for details). This vector field, in turn, induces a flow ψ^t on Γ_bound with respect to which one may define the time-translated test function:[Geometrically speaking, φ_t is the push-forward of φ with respect to ψ^t.]

    φ_t(x,p) := φ(ψ^-t(x,p)).

It is important to mention that in the rotating case, the Killing vector field ∂_t is spacelike and not timelike inside the ergoregion, meaning that (due to the dragging by rotation) stationary observers lying inside the ergoregion cannot follow the integral curves of ∂_t. In this case one might replace ∂_t with the linear combination X := ∂_t + Ω∂_φ, with the constant angular velocity Ω chosen such that X is timelike in the vicinity of the observer, and define ψ^t as the flow with respect to the complete lift X̂ of X. This provides a more sensible definition for the time-translated test function φ, if φ has its support in a region lying close to the event horizon which intersects the ergoregion. We stress that the theorem below holds for both (and probably more general) cases, the flow of ψ^t being defined with respect to the complete lift of either ∂_t or X.[As one can easily verify, the effect of replacing ∂_t by X is equivalent to replacing the frequency ω^1 by ω^1 + Ω  in the proof of Theorem <ref>.]

For the statement of the following theorem, the area function, defined by 

    A_m(E,L) := ∮ p_r dr = 2M m[ (1-ε^2)_0 + ε_3 ],
         L^2 > L_ms^2,   E_min(L) < E < E_max(L),

plays an important role. Note that this function determines the action variable J_2 = A_m/2π, and the frequencies ω^1 and ω^2 defined in Eq. (<ref>) are determined by the gradient of A_m as follows:

    ω^1 =  . ∂ A_m/∂ L/ ∂ A_m/∂ E.,
        ω^2 = -. 2π/ ∂ A_m/∂ E. .

After these remarks, we are ready to formulate the main result of this article:



Let F∈ L^1(S^1× S^1× (-∞,∞)× (-∞,∞)× (0,∞)) be a Lebesgue-integrable function which, according to Eq. (<ref>), determines a solution f(x,p) on Γ_bound of the collisionless Boltzmann equation on the equatorial plane of a Kerr black hole background. Let F and f be the corresponding distribution functions obtained by averaging over the angle variables. Let φ∈ C_0^∞(Γ_bound) be a test function, and denote by φ_t its time-translation as defined in Eq. (<ref>).

Suppose further that on the support of φ the following non-degeneracy condition holds:

    (D^2A_m(E,L)) ≠ 0,

where D^2 A_m(E,L) denotes the Hessian matrix of the area function (<ref>).

Then,

    lim_t→∞ N[φ_t] = ∫_Γ_boundf(x,p) φ(x,p)_Γ.



Remarks:


  * The validity of the determinant condition (<ref>) will be analyzed towards the end of this section, following the proof of the theorem. As will be verified numerically, it is satisfied everywhere except for points lying on a certain curve in (E,L)-space, see Fig. <ref>. However, we stress that the theorem does not require condition (<ref>) to hold everywhere; it is sufficient to hold for points lying in the support of φ, that is, in the vicinity where the observer performs the measurement.


  * Although the above formulation of the theorem requires the test function φ to be smooth, this assumption can be relaxed considerably. For example, it is possible to consider certain observables  with test functions of the form

    φ(x,p) = δ(x - x_obs)√(-(g^μν(x)))ψ(p),

corresponding to spacetime observables at some given event x_obs, as the ones considered in the previous section, see <cit.>.


  * The quantity n_∞ in Figs. <ref>, <ref> and <ref> has been calculated by replacing F by its average F in Eq. (<ref>). 


  * Note that the hypothesis of the theorem require φ to be supported on Γ_bound, which excludes circular orbits since the latter belong to the boundary of Γ_bound. However, the theorem can likely be generalized to the more realistic case in which the support of φ includes circular orbits (see also the remarks at the end of this section).




Theorem <ref> In a first step, we rewrite the integral defining N[φ] in terms of the action-angle variables (Q^α,J_α) as follows:

    N[φ] = ∫ F(Q^1 - ω^1 Q^0, Q^2 - ω^2 Q^0, J_0,J_1,J_2) 
    Φ(Q^0,Q^1,Q^2,J_0,J_1,J_2) dQ^0 dQ^1 dQ^2 dJ_0 dJ_1 dJ_2,

where Φ(Q^0,Q^1,Q^2,J_0,J_1,J_2) = φ(x,p) is the representation of the test function φ in terms of (Q^α,J_α) and where we have used the fact that the transformation (x,p)↦ (Q^α,J_α) is canonical. With respect to these variables, the flow associated with the vector field k̂ = -∂/∂ Q^0 amounts to a translation of Q^0 by -t, keeping the other variables fixed, such that φ_t(x,p) = Φ(Q^0 + t,Q^1,Q^2,J_0,J_1,J_2). By means of the variable substitution Q^0 = Θ^0 - t, Q^A = Θ^A + ω^AΘ^0, A = 1,2, one obtains

    N[φ_t] = ∫_0^∞ dm ∫_0^2π dΘ^1 ∫_0^2π dΘ^2 
    ∫_L^2 > L_ms^2 dL ∫_E_min(L)^E_max(L) dE  
    F_m(Θ^1 + ω^1 t, Θ^2 + ω^2 t, E,L) 
    Φ_m(Θ^1,Θ^2, E,L),

where we have defined F_m(Θ^1,Θ^2,E,L) := F(Θ^1,Θ^2,E,L,A_m(E,L)/2π) and

    Φ_m(Θ^1,Θ^2,E,L) 
     := 1/2π∂ A_m/∂ m(E,L)
    ∫_-∞^∞Φ( Θ^0,Θ^1 + ω^1Θ^0,Θ^2 + ω^2Θ^0,E,L,A_m(E,L)/2π)  dΘ^0.

In a next step, we perform a variable substitution in order to replace the integral over (E,L) by an integral over  (ω^1,ω^2). To this purpose, we denote by W_m: (E,L)↦ (ω^1,ω^2) the map that defines (for fixed m) the angular frequencies (ω^1,ω^2) in terms of the constants of motion (E,L). Its Jacobian determinant is related to the determinant of the Hessian of the area function as follows:

    (D W_m) = -(ω^2)^3/(2π)^2 (D^2A_m(E,L)).

Therefore, due to condition (<ref>) the map W_m is (at least) locally invertible on the support of φ. If not globally invertible, we may cover its domain with open subsets U_i on which W_m is invertible. Since Φ is compactly supported, only a finite number of these subsets are required. On each of these U_i's we define

    G_m,i(Θ^1,Θ^2,ω^1,ω^2) := {[ F_m(Θ^1,Θ^2,E,L))/| DW_m(E,L)| ,             (ω^1,ω^2)∈ W_m(U_i),;                               0,                       otherwise. ].

Assuming first that Φ_m is supported in one of these sets U_i, we have

    N[φ_t] =  ∫_0^∞ dm ∫_0^2π dΘ^1 ∫_0^2π dΘ^2 
    ∫_^2 dω^1 dω^2 g_m,i(t,Θ^1,Θ^2,ω^1,ω^2)
     Φ_m(Θ^1,Θ^2,E,L),

where

    g_m,i(t,Θ^1,Θ^2,ω^1,ω^2) := G_m,i(Θ^1 + ω^1 t,Θ^2 + ω^2 t,ω^1,ω^2).


After these initial steps, we encounter ourselves exactly in the same situation as a collisionless gas in a periodic box, and the mixing property is easily revealed by means of Fourier transformation (see section 3 in <cit.>): for this, define

    Ĝ_m,i(k_1,k_2,η_1, η_2) := 1/(2π)^2∫_-∞^∞∫_-∞^∞∫_0^2π∫_0^2πG_m,i(Θ ^1 , Θ ^2,ω ^1, ω ^2)
     e^-i(k_1 Θ ^1 + k_2 Θ ^2) e^- i(η_1ω^1 + η_2ω^2)  
      dΘ ^1 d Θ ^2 dω ^1 dω ^2 ,

where (k_1, k_2 ) ∈^2 and (η_1, η_2) ∈^2. Then,

    ĝ_m,i(t,k_1, k_2,η_1, η _2) = Ĝ_m,i(k_1,k_2,η_1 - k_1 t, η _2 - k_2 t),
 
hence the Fourier transform converts the rotations of the angle variables (Θ^1,Θ^2)↦ (Θ^1 + ω^1 t,Θ^2 + ω^2 t) into a translation (η_1,η_2)↦ (η_1 - k_1 t,η_2 - k_2 t) of the frequencies associated with the angular frequencies. According to the Riemann-Lebesgue lemma ĝ_m,i(t,k_1, k_2,η _1, η _2) converges pointwise to 0 for all fixed (k_1, k_2,η _1, η _2) with (k_1,k_2) ≠ (0,0). Consequently,

    lim_t→∞ĝ_m,i(t,k_1, k_2, η _1, η _2) = δ_k_1 0δ_k_2 0Ĝ_m,i(0,0,η_1, η _2) =: ĝ_m,i^∞(k_1, k_2,η _1, η_2),

for all (k_1, k_2,η _1, η _2)∈^2×^2 with the convergence rate determined by the smoothness of G_m,i(Θ ^1, Θ ^2,ω ^1, ω ^2) with respect to ω ^1, ω ^2. Here, ĝ_m,i^∞(k_1, k_2,η _1, η_2) is the Fourier transform of the function

    g^∞_m,i(Θ ^1, Θ ^2,ω ^1, ω^2) = 1/(2π)^2∫_0^2π∫_0^2π G_m,i(Θ ^1, Θ ^2,ω^1, ω^2) dΘ ^1 dΘ^2 
     = G_m,i(Θ ^1, Θ ^2,ω^1, ω^2),

which is the average of G_m,i over the angle variables. Since Φ_m is smooth and has compact support, it follows that N[φ_t] converges for t→∞ to the same expression as the one in the right-hand side of in Eq. (<ref>) with g_m,i(t,Θ^1,Θ^2,ω^1,ω^2) replaced with g^∞_m,i(Θ ^1, Θ ^2,ω ^1, ω^2). This proves the theorem for the case in which the support of Φ_m lies within one of the subsets U_i.

For the general case, a partition of unity can be used to write Φ_m as a finite sum of functions each of which is supported in only one of the subsets U_i.



Now that the theorem has been proven, we discuss the validity of the non-degeneracy condition (<ref>). Based on the explicit representation in Eq. (<ref>), we have computed the determinant of the Hessian of the area function. As in the previous section, it is convenient to describe the result in terms of the generalized Keplerian variables (P,e) instead of (E,L). Interestingly, it turns out the determinant condition is satisfied everywhere except for points lying on the blue solid curves shown in Figs. <ref> for different values of the rotational parameter.



To provide analytic support for the results shown in these figures we consider the particular case of quasi-circular orbits on a Schwarzschild background, for which α = 0 and e≪ 1. In this case, the expansion of the first derivatives of the area function in terms of the eccentricity e yield

    A_E    :=   ∂ A_m/∂ E = M2π P^2/√(P-6)[ 1 + 3/42P^3 - 32P^2 + 165P - 266/(P-2)(P-6)^2 e^2 +  O(e^4) ],
    
    
    A_L    :=   ∂ A_m/∂ L = -2π√(P)/√(P-6)[ 1 + 3/41/(P-6)^2 e^2 +  O(e^4) ],

giving

    ( [ ∂ A_E/∂ P ∂ A_E/∂ e; ∂ A_L/∂ P ∂ A_L/∂ e ])
     = M D(P) e +  O(e^3),    
     D(P) := -9π^2 P^3/2 (4P^2 - 39P + 86)/(P-2)(P-6)^3,     P > 6.

The function D is positive when P is slightly larger than 6 and negative for large P, and it has a single root at P = P_* := (39 + √(145))/8 ≃ 6.38, which corresponds to the limit point of the blue solid curve as e→ 0 in the middle panel of Fig. <ref>. From Eqs. (<ref>,<ref>) one can also conclude that (A_E,A_L) is a function of (P,e^2) which is locally invertible for small e^2 and P away from P_*.

We end this section by observing that the distribution functions considered in the previous section in cases A, E and F have their main support away from the blue solid curve in (P,e)-space where the determinant condition is violated, while in the remaining cases B, C and D the main support of the distribution function intersects this curve (compare the values given in Table <ref> with the middle panel of Fig. <ref>). Although the determinant condition is violated in the latter cases, the plots in Figs. <ref> and <ref> suggest that the particle density still converges for t→∞, indicating that mixing is still sufficiently strong for the relaxation process to take place.






§ CONCLUSIONS




In this work, we have shown that a relativistic, collisionless kinetic gas propagating in the equatorial plane of a Kerr black hole spacetime settles down to a stationary, axisymmetric configuration. As we have demonstrated, this relaxation process is due to phase space mixing, an effect that plays a prominent role in a wide range of fields in physics including galactic dynamics and plasma physics. Here, we have exhibited the relevance of the mixing phenomenon for the dynamical behavior of spacetime observables within the fully general relativistic setting of a kinetic gas which is trapped in the strong field gravitational field of a black hole.

The main implication of this work is that the one-particle distribution function f describing the state of the gas, which in general is a function of the five coordinates (Q^1,Q^2,J_0,J_1,J_2) parametrizing the relativistic phase space Γ_bound corresponding to bound trajectories, can be replaced by a much simpler distribution function f which is a function depending only of the constants of motion (J_0,J_1,J_2) and which can be computed by averaging the initial distribution function over the angle variables (Q^1,Q^2). Indeed, our main theorem in the previous section shows that, provided the determinant condition (<ref>) is satisfied on the support of φ, the integral over f times any test function φ converges in time to the integral over f times the same test function φ. At the physical level, these integrals have the interpretation of macroscopic observables, where different choices for the test function correspond to different physical quantities being measured.

The determinant condition (<ref>) means that the map (E,L)↦ (ω^1,ω^2) which defines the angular frequencies in terms of the constants of motion E and L is locally invertible, a condition that is well-known in perturbation theory of integrable Hamiltonian systems, see for example Section X.51 in <cit.>. For geodesic motion in the equatorial plane of a Kerr black hole we have shown that the determinant condition holds everywhere with the exception of points in the (E,L)-plane lying on a curve (a zero-measure set) which connects a particular circular orbit to innermost stable orbits at large eccentricities. As the plots in Section <ref> indicate, mixing still occurs in the vicinity of these exceptional points, suggesting that our theorem also holds under weaker assumptions. It should be interesting to relate the behavior of the gas in the vicinity of the exceptional points to Tremaine's analysis of stable singularities or catastrophes in galaxies <cit.> and the implications on the time scale in which the mixing occurs. In any case, for the specific examples we have analyzed in Section <ref>, the damping of the oscillations in the particle density is rather fast, with relative amplitude lying below 0.001 after a few thousand light-crossing times corresponding to the gravitational radius of the black hole.

A further interesting problem consists in analyzing the effects of the self-gravity of the gas configuration (which have been neglected in the present work) on the mixing property. The inclusion of the self-gravity implies that the Kerr metric acquires correction terms due to the non-vanishing stress-energy tensor associated with the gas, which in turn leads to a perturbed Hamiltonian flow describing the geodesic motion for the gas particles. Based on general arguments from Kolmogorov-Arnold-Moser (KAM) theory, it has been argued <cit.> that such perturbations could trigger dynamical chaos in the vicinity of resonant orbits. Therefore, it should be particularly interesting to study the mixing phenomenon in the neighborhood of such orbits, and investigate whether or not the relaxation of the observables persists in the self-gravitating case. 

The generalization of the mixing property to thick disk configurations, where individual gas particles are not necessarily confined to the equatorial plane, will be given in future work <cit.>.








We wish to thank Francisco Astorga, Michael Kiessling, Alexander Komech, Ulises Nucamendi, and Thomas Zannias for fruitful and stimulating discussions. We also thank Michael Kiessling, and Thomas Zannias for comments on a previous version of this manuscript. This research was supported in part by CONACyT Grants No. 577742, by the CONACyT Network Project No. 294625 “Agujeros Negros y Ondas Gravitatorias", and by a CIC Grant to Universidad Michoacana. We also thank the Erwin Schrödinger Institute for Theoretical Physics, where this work was initiated, for hospitality. P.R. also thanks Chalmers University of Technology and the University of Vienna for hospitality, where part of this work was completed.






§ EXPLICIT EXPRESSIONS FOR THE FUNCTIONS _J




Using the definition of the function R(r) in Eq. (<ref>), the following relations between (ε,λ) and its roots (ξ_j) are obtained:

    1 - ε^2 = 2/ξ_012,    λ̂^2 = ξ_0 ξ_1 ξ_2/ξ_012,    λ̂^2 + 2αλ̂ε + α^2 = 2ξ_0 ξ_1 + ξ_0 ξ_2 + ξ_1 ξ_2/ξ_012,

where we have abbreviated ξ_012 := ξ_0 + ξ_1 + ξ_2 and introduced λ̂ := λ - αε. The functions  _j(χ) and corresponding constants _j := _j(π/2) in terms of which the angle variables Q^α are expressed are defined as follows:

    _0(χ)    :=    -C/2{(ξ_0 ξ_012 - ξ_1 ξ_2)(χ,k)
     + ξ_1(ξ_2 - ξ_0)(χ,k) .
          . + (ξ_1 - ξ_0) ξ_012Π(χ,b^2,k) - 1/2ξ_1(ξ_2 - ξ_1)
    √(1 - k^2sin^2χ)/1 - b^2sin^2χsin(2χ) },
    
    _1(χ)    :=    -C{( λ̂ + αεξ_0^2 - αλ̂/(ξ_0 - ξ_+)(ξ_0 - ξ_-))(χ,k)
    .
          . - αξ_1 - ξ_0/ξ_+ - ξ_-[ εξ_+^2 - αλ̂/(ξ_0 - ξ_+)(ξ_1 - ξ_+)Π(χ,b_+^2,k) - (+↔ -) ] },
    
    _2(χ)    :=    2C{ξ_0(εξ_0^2 - αλ̂)/(ξ_0 - ξ_+)(ξ_0 - ξ_-)(χ,k) + ε(ξ_1 - ξ_0)Π(χ,b^2,k) .
          . - ξ_1 - ξ_0/ξ_+ - ξ_-[ ξ_+(εξ_+^2 - αλ̂)/(ξ_0 - ξ_+)(ξ_1 - ξ_+)Π(χ,b_+^2,k) - (+↔ -) ] },

and

    _3(χ)    :=   _2(χ) + λ/ε_1(χ)
    
       =   2C/√(ξ_012(ξ_012-2)){ξ_0^2(χ,k) + (ξ_1 - ξ_0)(ξ_012 - 2)Π(χ,b^2,k) .
          . - ξ_1 - ξ_0/ξ_+ - ξ_-[ ξ_+(ξ_2 - ξ_+)Π(χ,b_+^2,k) - (+↔ -) ] },

with b := √((ξ_2 - ξ_1)/(ξ_2 - ξ_0)), k := √(ξ_0/ξ_1) b, ξ_± := 1 ±√(1 - α^2), b_± := √(ξ_0 - ξ_±/ξ_1 - ξ_±) b and C := √(2ξ_012/ξ_1(ξ_2-ξ_0)), and where the functions (χ,k), (χ,k) and Π(χ,b^2,k) denote standard elliptic integrals as defined, for instance, in Ref. <cit.>.

In terms of these functions one finds

    S = M m[ -εt/M + λφ +  (1-ε^2)_0(χ) + ε_3(χ) ],

and

    ∂ S/∂ m = M_0(χ),    ∂ S/∂ L = φ + _1(χ),    ∂ S/∂ E = -t + M_2(χ) - E_0(χ),

from which the expressions (<ref>) are easily derived.




unsrt


10

dL62
D. Lynden-Bell.
The stability and vibrations of a gas of stars.
Royal Astronomical Society, 124:279–296, 1962.

dL67
D. Lynden-Bell.
Statistical mechanics of violent relaxation in stellar systems.
Monthly Notices Roy Astronom. Soc., 136:101–121, 1967.

sTmHdL86
S. Tremaine, M. Hénon, and D. Lynden-Bell.
H-functions and mixing in violent relaxation.
Monthly Notices Roy Astronom. Soc., 219:285–297, 1986.

sT99
S. Tremaine.
The geometry of phase mixing.
Monthly Notices Roy Astronom. Soc., 307:877–883, 1999.

BinneyTremaine-Book
J. Binney and S. Tremaine.
Galactic Dynamics (Second Edition).
Princeton University Press, Princeton, New Jersey, 2008.

gCrSmFbGpKpA14
G. N. Candlish, R. Smith, Fellhauer, B. K. Gibson, P. Kroupa, and P. Assmann.
Phase mixing due to the galactic potential: steps in the position and
  velocity distributions of popped star clusters.
Monthly Notices Roy Astronom. Soc., 437:3702–3717, 2014.

cMcV11
C. Mouhot and C. Villani.
On Landau damping.
Acta Math., 207:29–201, 2011.

bY15
B. Young.
On linear Landau damping for relativistic plasmas via Gevrey
  regularity.
J. Diff. Eqns., 259:3233 – 3273, 2015.

bY16
B. Young.
Landau damping in relativistic plasmas.
J. Math. Phys., 57:021502, 2016.

rMeT17
R. Mathew and E. Tiesinga.
Phase-space mixing in dynamically unstable, integrable few-mode
  quantum systems.
Phys. Rev. A, 96:013604, 2017.

tDaKeKyS02
T.V. Dudnikova, A. I. Komech, E. A. Kopylova, and Y. M. Suhov.
On convergence to equilibrium distribution, I. the KleinGordon
  equation with mixing.
Comm. Math. Phys., 225:1–32, 2002.

tDaKnRyS02
T.V. Dudnikova, A. I. Komech, N.E. Ratanov, and Y. M. Suhov.
On convergence to equilibrium distribution, II. the wave equation
  in odd dimensions, with mixing.
J. Stat. Phys., 108:1219–1253, 2002.

pRoS18c
P. Rioseco and O. Sarbach.
Dynamics of collisionless thick discs around a Kerr black hole.
in preparation, 2018.

wS02
W. Schmidt.
Celestial mechanics in Kerr space-time.
Class. Quantum Grav., 19:2743–2764, 2002.

tHeF08
T. Hinderer and E.E. Flanagan.
Two timescale analysis of extreme mass ratio inspirals in Kerr. I.
  Orbital Motion.
Phys. Rev. D, 78:064028, 2008.

rFwH09
R. Fujita and W. Hikida.
Analytical solutions of bound timelike geodesic orbits in Kerr
  spacetime.
Class. Quant. Grav., 26:135002, 2009.

jBwPsT72
J. Bardeen, W. Press, and S. Teukolsky.
Rotating black holes: Locally nonrotating frames, energy extraction,
  and scalar synchrotron radiation.
Astrophysical Journal, 178:347–370, 1972.

Arnold-Book
V.I. Arnold.
Mathematical Methods of Classical Mechanics.
Springer-Verlag, New York, 1989.

oStZ14b
O. Sarbach and T. Zannias.
The geometry of the tangent bundle and the relativistic kinetic
  theory of gases.
Class. Quantum Grav., 31:085013, 2014.

pRoS16
P. Rioseco and O. Sarbach.
Accretion of a relativistic, collisionless kinetic gas into a
  Schwarzschild black hole.
Class. Quantum Grav., 34(9):095007, 2017.

jE71
J. Ehlers.
General relativity and kinetic theory.
In R.K. Sachs, editor, General Relativity and Cosmology, pages
  1–70, New York, 1971. Acedemic.

jBmGtH15
J. Brink, M. Geyer, and T. Hinderer.
Astrophysics of resonant orbits in the Kerr metric.
Phys. Rev. D, 91(8):083001, 2015.

pRoS18b
P. Rioseco and O. Sarbach.
Phase space mixing in external gravitational potentials.
in preparation, 2018.

jBmGtH13
J. Brink, M. Geyer, and T. Hinderer.
Orbital resonances around black holes.
Phys. Rev. Lett., 114:081102, 2015.

DLMF
Digital library of mathematical functions.
<http://dlmf.nist.gov/>.






