On the spatially homogeneous and isotropic Einstein-Vlasov-Fokker-Planck system with cosmological scalar field

The Einstein-Vlasov-Fokker-Planck system describes the kinetic diffusion dynamics of self-gravitating particles within the Einstein theory of general relativity. We study the Cauchy problem for spatially homogeneous and isotropic solutions and prove the existence of both global solutions and solutions that blow-up in finite time depending on the size of certain functions of the initial data. We also derive information on the large-time behavior of global solutions and toward the singularity for solutions which blow-up in fine time. Our results entail the existence of a phase of decelerated expansion followed by a phase of accelerated expansion, in accordance with the physical expectations in cosmology.


Introduction
The purpose of this paper is to study spatially homogeneous and isotropic solutions of the Einstein-Vlasov-Fokker-Planck system. The model describes the kinetic diffusion dynamics of self-gravitating particles within the Einstein theory of general relativity. It is assumed that diffusion takes place in a cosmological scalar field, which can be identified with the dark energy source responsible for the phase of accelerated expansion of the Universe.
When relativistic effects are neglected the motion of self-gravitating kinetic particles undergoing diffusion is described by the frictionless Vlasov-Poisson-Fokker-Planck system in the gravitational case, which is given by f (t, x, p) dp. (1b) Here f (t, x, p) is the phase-space density of a system of unit mass particles, U (t, x) is the Newtonian gravitational potential generated by the particle system, and σ > 0 is the diffusion constant; the remaining physical constants have been set to one. The mathematical properties of the system (1) have been extensively studied in the literature. In particular, it is known that the Vlasov-Poisson-Fokker-Planck system admits a global unique solution for given initial data [4,11], whose asymptotic behavior for large times has been studied in [6]. The proof of both these results makes use of the explicit form of the fundamental solution to the linear Fokker-Planck equation. Similar results have also been obtained for the related Vlasov-Maxwell-Fokker-Planck system [12,13]. A relativistic generalization of the Vlasov-Poisson-Fokker-Planck system in the gravitational case has been introduced recently in [5] and the purpose of this paper is to initiate its mathematical study. When dealing with the relativistic model one is faced with many new difficulties, including the hyperbolic and nonlinear character of the Einstein field equations (compared to the linearity of the Poisson equation in (1)), the non-uniform ellipticity of the diffusion operator, the time-dependence of the diffusion matrix, and the absence of an explicit formula for the fundamental solution to the linear relativistic Fokker-Planck equation. The further notorious complexity of the Einstein equations suggests that one begins by studying solutions with symmetries. In this paper we consider solutions which are spatially homogeneous and isotropic. Under these symmetry assumptions the Einstein equations reduce to a system of nonlinear ordinary differential equations.
In the absence of diffusion our model reduces to the Einstein-Vlasov system with cosmological constant. The regularity and asymptotic behavior of spatially homogeneous solutions to the latter system have been studied in [8]. The generalization of these results to spatially inhomogeneous solutions (with surface symmetry) is given in [3]. When the cosmological constant is set to zero one obtains the Einstein-Vlasov system. Some results available for this system are summarized in [1]. It has been shown recently in [7,9] that the Cauchy problem for the Einstein-Vlasov system is globally well-posed for small data with no symmetry restrictions.
In [2] we have studied a similar problem when the Einstein equation is replaced by a nonlinear wave equation for a scalar field. It turns out that for spatially homogeneous and isotropic solutions, the Fokker-Planck equation considered in [2] is the same as the one derived in the present paper and thus the analysis of the matter equation does not pose any new difficulty. Contrastingly, the Einstein equations behave quite differently than the field equation considered in [2] The remainder of the paper is organized as follows. In the next section we derive the Einstein-Vlasov-Fokker-Planck system for spatially homogeneous and isotropic solutions. In Section 3 we study the Cauchy problem and prove that, depending on the size of certain functions of the the initial data, regular solutions either exist globally or blow-up in finite time. In Section 4 we derive information on the asymptotic behavior as t → ∞ of global solutions and toward the singularity for solutions which blow-up in finite time. We also show that for an open set of initial data which describe an initially decelerating expanding Universe, global solutions will eventually give rise to a phase of accelerated expansion, in agreement with the expectations in Cosmology [14].

Derivation of the model
We begin with a short description of the general relativistic kinetic theory of diffusion, see [5] for more details. In the following, Greek indices run from 0 to 3, while Latin indices run from 1 to 3. Moreover all indices are raised and lowered with the matrix g µν and its inverse g µν := (g −1 ) µν , and the Einstein summation rule is applied (e.g., g µν g να = δ µ α ), unless otherwise stated. Let (M, g) be a spacetime. Let x µ denote local coordinates on M , with t = x 0 being timelike, and let (x µ , p ν ) be the induced canonical local coordinates on the tangent bundle of the spacetime. The mass-shell for particles with unit mass is the 7-dimensional submanifold of the tangent bundle obtained by imposing g µν p µ p ν = −1 and p 0 > 0, so that (p 0 , p 1 , p 2 , p 3 ) can be interpreted as the components of the four-momentum of a unit mass particle. The mass-shell condition is used to express p 0 in terms of the spatial components of the four-momentum, namely It follows that (x µ , p i ) define a local system of coordinates on the mass-shell. The particle distribution function f is defined on the mass-shell and is therefore a function of the coordinates (t, x i , p j ). The Fokker-Planck equation for f is the PDE where L is the Liouville (or Vlasov) operator, Γ α µν denote the Christoffel symbols of the metric g, σ is the (positive) diffusion constant, and D p is the diffusion operator, which is the Laplace-Beltrami operator associated to the Riemannian metric h = h ij dp i dp j , h ij = g ij + g 00 We have where (h −1 ) ij is the inverse matrix of h ij (i.e., (h −1 ) ij h jk = δ i k ). From the particle distribution f we can construct the particle current density and the energy-momentum tensor by where dp = dp 1 ∧ dp 2 ∧ dp 3 and the integration is over the fibers of the mass-shell. It can be shown that the tensors T µν and J µ verify To close the system we require the metric g to solve the Einstein equations with cosmological scalar field φ, which, in units 8πG = c = 1, is given by By (6) and the Bianchi identity ∇ µ (R µν − 1 2 g µν R) = 0, we find that the cosmological scalar field satisfies the equation In particular, φ satisfies the wave equation ✷φ = ∇ µ ∇ µ φ = 3σ∇ µ J µ = 0 and where for the inequality we have used the fact that J µ is timelike. Hence, the cosmological scalar field is decreasing along the matter flow, which can be interpreted as energy being transferred from the scalar field to the particles by diffusion. Next, we specialize to spatially homogeneous and isotropic solutions. Under these symmetry assumptions, the spacetime metric can be written as where K(r) = 1 + k 4 r 2 , r = |x|, k ∈ {−1, 0, 1}.
Here x = (x 1 , x 2 , x 3 ) is a system of spatial isotropic coordinates on the hypersurfaces t = const., with t denoting the proper time along the normal geodesics, and | · | denoting the standard Euclidean norm. The hypersurfaces of constant proper time have zero/positive/negative constant curvature according to the values k = 0, +1, −1 of the curvature parameter k. Equation (2) gives and p 0 = −p 0 . The non-zero Christoffel symbols Γ i µν of the metric (9) are given by and by their symmetric symbols on the lower indexes; the index i is not summed in the previous equations. It follows that the Liouville operator in the left hand side of (3) is given by with w · z denoting the standard Euclidean scalar product of the vectors w, z ∈ R 3 . The metric (4) now takes the form It follows that Hence, the diffusion term in the right hand side of (3) ultimately takes the form Definition 2.1. A particle distribution f is said to be spatially homogeneous and isotropic if there exists a function F : It can be shown that this definition is equivalent to require that f is invariant by the six-dimensional group of isometries of the metric (9), see [10]. It is now convenient to define the function F : in terms of which the Fokker-Planck equation for spatially homogeneous and isotropic distribution functions takes the final form where the diffusion matrix D is In terms of the function F , and recalling (12), the energy-momentum tensor and current density (5) read and T 0i = T ij = 0, for i = j. Thus, equation (8) for the cosmological scalar field where is the total number of particles, which is conserved along solutions of (13).
Finally the non-zero components of the Einstein tensor G µν = R µν − 1 2 g µν R are It follows that the Einstein equations (7) are

By introducing the Hubble function
as well as the energy density ρ(t) and the pressure P(t) by we can rewrite the Einstein equations in the form We also have the following auxiliary equations, obtained by combining (20): Note also the estimates which follow straightforwardly by the definitions (18) and (19). The deceleration parameter Q is defined as The solution is said to be expanding with acceleration if Q < 0 and H > 0 and expanding with deceleration if Q > 0 and H > 0. From (21), this occurs at time t if and only if H(t) > 0 and is negative, respectively positive.

Existence of regular solutions
In this section we study the existence of regular solutions to the initial value problem for the spatially homogeneous and isotropic Einstein-Vlasov-Fokker-Planck system with cosmological scalar field. The system, derived in the previous section, is given by where ρ(t), P(t) are given by (19) and k is either 0, +1, or −1. Initial data for the system (26) consist of a quadruple We assume that F 0 is not identically zero and belongs to the space Moreover, the initial data are assumed to satisfy (26e) at time t = 0, namely Given T > 0, a quadruple (F (t, v), a(t), H(t), φ(t)) will be referred to as a regular solution of the system (26) on the interval [0, T ) and with initial data Note that for regular solutions the functions ρ, P are continuous and hence the Einstein equations are satisfied in the pointwise, classical sense. However the Fokker-Planck equation need only be verified in the weak sense.
Proving existence and uniqueness of local regular solutions is a simple generalization of the argument presented in [2]. The main tool is the following set of a priori estimates on the solution of the Fokker-Planck equation. The following estimates on F hold for all t ∈ [0, T max ): where α = α(σ) > 0, with α(0) = 0, and C > 0 are constants (depending on γ) and (y) + denotes the positive part of y.
Proof. The proof is very similar to the proof of Prop. 2.2 in [2], hence we limit ourself to formally derive the estimate in (ii). We compute In the last integral we use the brief calculation and Grönwall's Lemma concludes the proof of (ii).
The following lemma provides a characterization of ρ(t) that will be useful in subsequent results.
In order to determine whether a regular solution exists globally or blows-up in finite time, we shall often apply the following simple result within subsequent sections. Lemma 3.3. Let t 1 ∈ R be given and assume there is t 2 > t 1 with I ∈ C 1 (t 1 , t 2 ) satisfying for all t ∈ (t 1 , t 2 ).
(a) If t 2 ≥ t 1 − 1 I 1 , then t 2 = t 1 − 1 I 1 and lim t→t − 2 I(t) = −∞ with the estimate Proof. Sinceİ(t) ≤ 0 on (t 1 , t 2 ) and I 1 < 0, it follows that I(t) < 0 on the same interval and we divide by −I(t) 2 so that the differential inequality becomes for all t ∈ (t 1 , t 2 ). Integrating over [t 1 , τ ] for τ ∈ (t 1 , t 2 ) yields and for τ ∈ t 1 , t 1 − 1 I 1 , this can be inverted to find Taking the limit as τ → t 1 − 1 implies I(τ ) → −∞ at this time. Finally, since the differential inequality cannot be satisfied after the blow-up time t 1 − 1 I 1 , it follows that t 2 = t 1 − 1 I 1 . The resulting estimate then follows by making this replacement in the above inequality.
In the second case, we first note that the limit condition implies lim t→t − 2 I(t) −1 = 0. Integrating (28) over [τ, t 2 ) for τ ∈ (t 1 , t 2 ), we find and thus the lower bound

Characterization of the maximal time of existence
Next, we focus on proving a series of criteria for the existence of global regular solutions, or their finite-time blowup. We begin by showing that as long as a(t) remains bounded away from zero, the solution remains regular.
As H =ȧ/a we also have a ∈ W 1,∞ ([0, T max )) and the proof is complete.
Theorem 3.5. Let T max > 0 be the maximal time of existence of a regular solution.
Then, for k = 0, −1 the following are equivalent: For k = 1, we can only prove the weaker statements (b)⇒(c): As H(t) is strictly positive for t ∈ [0, T max ), then a(t) is increasing and so a(t) ≥ a 0 for all t ∈ [0, T max ). Statement (c) then follows immediately by Lemma 3.4.
Finally, we complete the proof by showing that either positive lower bound on φ in the case k = 1 implies global existence. First, assume As before, since φ is decreasing, we have φ(t) > 4 N 2 for all t ∈ [0, T max ). Using the lower bound ρ(t) ≥ N/a(t) 3 in (26e) we obtain Defining the function for all x > 0, we find Therefore, g ′ (x) = 0 only at x = N 2 , and g is minimized at this point. Thus, This implies that H does not change sign and thus H(t) > 0 on [0, T max ). Hence, a(t) > 0, which implies a(t) ≥ a 0 for all t ∈ [0, T max ), and T max = ∞ follows from Hence, a lower bound on ρ follows, namely Thus, on the same time interval, we use (26e) to find Defining the function for x > 0, we find its minimum value occurs when x 2 = 2 3 ρ 0 a 4 0 , and at this value Using the above assumption on the strictly decreasing function φ(t), we see that H(T * ) remains positive, and as in the previous argument this implies T * = T max = ∞.

Initial data for global existence and blowup
Our next goal is to prove that for each k = 0, ±1 there exist conditions on the initial data for which the solution is global and conditions which imply finite time blow-up. There are many ways to express these conditions; one is to write them in the form φ 0 < · · · ⇒ blow-up and φ 0 > · · · ⇒ global existence, where the right hand side of the inequality depends on N, H 0 , a 0 and the given constants k and σ (but not on φ 0 ). When expressed in this form it is straightforward to show that the conditions on the initial data derived in the present section are compatible with the constraint equation (27). We begin with the blow-up results.
Theorem 3.6. Let k = 0, −1. Then for the regular solution blows-up in finite time, i.e., T max < ∞. If k = 1, then the solution blows-up in finite time when where Erfc(z) denotes the complementary error-function.
In the case k = 1, the bound on H is much weaker. Hence, let us suppose that (33) holds and assume T max = +∞. We first notice that by (26d) and the bound Integrating (26d) we obtain as t → ∞. Thus, under condition (33) we have lim t→+∞ φ(t) < 0, and by Theorem 3.5 we reach a contradiction.
With this, we turn to initial data which launch a global-in-time solution.
Using (26d) we get Hence, under condition (34), if T * < T max we would have φ(T * ) > 0, which contradicts the maximality of T * . It follows that T * = T max = +∞. The lower bound on the limit of φ(t) follows by letting t → ∞ in (35).
Notice that the last two results suggest a particular quantity involving initial data that can predict the finite or infinite lifespan of the smooth solution in the cases of k = −1, 0. More specifically, if φ 0 H 0 a 3 0 σN is large enough then the solution must be global, while if this quantity is made too small, then solutions must blow-up in finite time. As before, the situation for k = 1 is more complicated, and the next result establishes conditions on initial data that guarantee global existence in this case.

Now, the function
is increasing for all t ≥ 0 and the minimum for nonnegative times occurs at t = 0. Therefore, assuming that (36) holds. Instead, if 3σ H 0 > 1 then t * > 0. Therefore, assuming that (37) holds. In either case,Ḣ(t) > 0 for all t ∈ (0, T * ), which implies that H(t) > H 0 on the same interval. Hence, the above lower bound on a(t) continues on [0, T * ), yielding a regular solution, and because H(T * ) > H 0 we find T * = ∞.
Our final result in this section is meant to unify the treatment for all three cases. and suppose that (38) holds. After some algebra, we see that this inequality implies Then, define and note that since β < 1, we have T * > 0. Now, integrating (26d) for t ∈ [0, T * ), we find Hence, (26e) implies so that H(t) does not change sign on [0, T * ) and H(t) ≥ γ. It follows that a(t) ≥ a 0 e γt for all t ∈ [0, T * ). Since β < 1, this then implies T * = T max . Finally, since γ > 0 we find inf It follows by Lemma 3.4 that T max = ∞, and the proof of the theorem is complete.

Asymptotic behavior of solutions
In this section we analyze the asymptotic behavior as t → ∞ of global solutions and as t → T − max of solutions which blow-up in finite time. We begin with the former solutions. Theorem 4.1. Consider a global regular solution, i.e., T max = ∞, and let φ ∞ := lim t→∞ φ(t) ≥ 0. Then for k = 0, −1 the following holds: (ii) If φ ∞ = 0, then there exist C 1 , C 2 > 0 such that for all t ≥ 0 4ρ 0 a 4 0 := φ m then the same conclusion as in (i) holds with φ ∞ replaced by φ * = φ ∞ − φ m .
Proof. We first assume k = 0, −1. If φ ∞ > 0, then, by (26e), for all t ≥ 0. Thus H remains strictly positive for all time with H(t) ≥ φ∞ 3 . It follows that for all t ≥ 0 Furthermore, due to (19) and using the positivity of H in Lemma 3.2 implies N a(t) −3 ≤ ρ(t) ≤ (ρ 0 a 3 0 + 3σN t)a(t) −3 which yields the claimed behavior of ρ once the bounds on a(t) are established. To arrive at a bound on φ, we integrate (26d) and use the lower bound on a to find for all t ≥ 0, which implies an upper bound on the behavior of φ. Using this in (26e) with the upper bound on ρ and the lower bound on a implies for all t ≥ 0 This further implies and, in conjunction with (39), this inequality implies the claimed behavior of H. Finally, this upper bound on H implies an analogous upper bound on a, so that Pairing this with (40) implies the bounds on a(t), and the upper bound on a(t) further implies the lower bound on φ(t) using (26d). Next, assume φ ∞ = 0. As in the previous case, we find by (26e) for all t ≥ 0, and thus H remains strictly positive for all time. As in the proof of Theorem 3.7, we use (26e), (21), and the bound P ≤ ρ/3 to finḋ for t ≥ 0. Integrating and using the positivity of H(t) we obtain and the lower bound follows for t ≥ 0. Furthermore, the positivity of H implies the same upper bound on ρ as in the previous case, and we again find for t ≥ 0 As before, this yields the claimed asymptotic behavior of ρ once the full behavior of a is obtained. Integrating (26d) and using (42), we find Next, we use (26e), the upper bounds on ρ and φ, and the lower bound on a so that and thus With this, an upper bound on a(t) follows, namely The last claim for k = 1 follows by the same argument used to prove (i) and the lower bounds on H attained for k = 1 within the proof of Theorem 3.5. For example, assuming φ ∞ ≥ 4 N 2 the argument proceeds by using instead of (39). Similarly, assuming φ ∞ ≥ 9 4ρ 0 a 4 0 the argument proceeds by using instead of (39).
We remark that sufficient conditions to have φ ∞ > 0 for k = 0, −1 and φ * > 0 for k = 1 follow by the lower bounds on lim t→∞ φ(t) derived in Theorems 3.7 and 3.9. Next, we provide some information on the behavior of solutions that blow-up in finite time as they approach the singularity. Moreover, the singularity is a curvature singularity, and there exist C 1 , C 2 > 0 such that for t sufficiently close to T max For the case k = 1, if lim Proof. By Theorem 3.5 the finite time blow-up of the solution in the cases k = 0, −1 implies H(t) ≤ 0 for some t ∈ (0, T max ) and φ(t) must eventually attain negative values. Without loss of generality, we can take T 1 ∈ (0, T max ) large enough so that φ(t), H(t) < 0 for t ∈ [T 1 , T max ), as φ is decreasing and φ(t) < 0 impliesḢ(t) < 0 by (21). With this, we see that a(t) is strictly decreasing on (T 1 , T max ), and by Lemma 3.4 we find both a(t) > 0 for all t ∈ [0, T max ) and inf t∈[0,Tmax) a(t) = 0.
Thus, it follows that lim t→T − max a(t) = 0. The limits of H and ρ in (43) follow from this behavior. Indeed, the limiting behavior of ρ follows directly from the lower bound (23). Also, because φ(t) < 0 for t ∈ [T 1 , T max ), we use (21) to arrive aṫ for t ∈ [T 1 , T max ). Then, using the definition of H, we computė on the same time interval. Using this with the lower bounds on a and φ, we have for t sufficiently close to T max .
Next, we establish bounds on these functions in the opposite directions. Since we know a → 0 as t → T − max and H(T 1 ), φ(T 1 ) < 0 with each of these functions decreasing on (T 1 , T max ), we can find T 2 ∈ [T 1 , T max ) such that H(t), φ(t) < 0 and a(t) < 1 for t ∈ [T 2 , T max ). Then, using (22), (23), and (26d) we finḋ ρ(t) ≤ −4H(t)ρ(t) + 3σN a(t) −3 ≤ −4H(t)ρ(t) + 3σN a(t) −4 for t ∈ [T 2 , T max ). Multiplying by a(t) 4 , we can rewrite the inequality as for t ∈ [T 2 , T max ). Using this estimate within (26e), we find for t ∈ [T 2 , T max ). Multiplying by a(t) 4 , this becomes |a(t)ȧ(t)| 2 ≤ C which, becauseȧ(t) < 0, further implies Integrating over [τ, T max ) and using the limiting behavior of a as t → T − max , we find With the upper bound on a established, the lower bound on ρ follows directly from (23) so that for t ∈ [T 2 , T max ). Additionally, the upper bound on φ follows exactly as before. In particular, integrating (26d) over [T 2 , t] we find for t sufficiently close to T max . Furthermore, this estimate implies the limiting behavior of φ within (43). Finally, the upper bound on H is obtained from (21) so thaṫ for t ∈ [T 2 , T max ). Integrating over [T 2 , τ ] then yields the result for τ sufficiently close to T max . For the case k = 1 we merely repeat this argument under the assumption that φ(t) < 0 for some t ∈ (0, T max ) since this gives rise to negative values of H(t) as in the proof of Theorem 3.5, as well as, the upper boundḢ(t) ≤ −H(t) 2 .
With these results at hand we can finally discuss the important question of the existence of a phase of accelerated expansion of the Universe in the future of a phase of decelerated expansion. There is overwhelming experimental evidence that the Universe is currently expanding with acceleration. Nonetheless, the standard cosmological models require the existence in the past of a phase of decelerated expansion, during which the structures visible today (galaxies, clusters, etc.) have formed. The following last result shows that the diffusion model studied in this paper is able to reproduce this physical behavior.