A non-relativistic Model of Plasma Physics Containing a Radiation Reaction Term

While a fully relativistic collisionless plasma is modeled by the Vlasov-Maxwell system a good approximation in the non-relativistic limit is given by the Vlasov-Poisson system. We modify the Vlasov-Poisson system so that damping due to the relativistic effect of radiation reaction is included. We prove the existence and uniqueness as well as the higher regularity of local classical solutions. These theorems also include the higher regularity of classical solutions of the Vlasov-Poisson system depending on the regularity of the initial datum.


Introduction
In studying interactions of charged matter and electromagnetic fields it is a quite common situation that the velocity of matter is small if compared to the velocity of light. In such situations an asymptotic expansion of the equations in the small parameter |v| c can help to reduce the complexity of the original problem, e.g. for numerical computations. Since in many relativistic settings the zeroth order contribution of this expansion gives the Newtonian limit such expansions are usually called post-Newtonian expansions.
In numerous applications of asymptotic analysis it is quite straightforward to establish the equations of the expansion and it is the main task to prove a relation between the solutions of the asymptotic expansion and the solution of the original equation. In contrast to that in interactions of charged matter and electrodynamic fields also the formulation of the approximation scheme becomes delicate if effects due to radiation damping shall be included. It is well-known that accelerated charged matter loses energy by radiation; an electromagnetic radiation field is generated which transports energy to null-infinity. This gives a damping effect on the motion of the matter known as radiation damping. There is a large amount of literature for single particle models concerning effective equations which incorporate radiation damping without giving a fully relativistic description, see [19] and the literature cited therein.
On the other hand the Vlasov-Maxwell system is a well-known fully relativistic description of a large ensemble of collisionless particles which interacts by means of the electromagnetic fields collectively generated by these particles. In a non-relativistic setting this situation is governed by the Vlasov-Poisson system. While the asymptotic limit relation between the Vlasov-Poisson and the Vlasov-Maxwell system is well understood, see [18,14], and also higher order relations between the so called Vlasov-Darwin and the Vlasov-Maxwell model are established, see [3], all these approximation models do not inhibit radiation damping.
It is the aim of this paper to establish and study a model in between the fully relativistic description and the non-relativistic description where effects due to radiation damping are included. In Section 2 we repeat some heuristics from [11,12] leading to two candidates of such systems, see equations (6) and (7). While the first of these systems is already studied in [11,12] we give an additional heuristic leading to (7). In the third section we shall reformulate (7) establishing the topic of this paper, the 'reduced Vlasov-Poisson system with radiation damping', see (rVPRD ε ).
In this situation the main criterion to decide wether or not a model is suitable, is the approximation property. In a second paper [2] it is proved that solutions of (rVPRD ε ) together with some lower order terms approximate solutions of the Vlasov-Maxwell system to a higher order in v/c than solutions of the Vlasov-Poisson or the Vlasov-Darwin model, see also [4] for an overview.
In Section 4 we present a proof of local existence, uniqueness and smoothness of solutions of (rVPRD ε ). This proof is an adaption of the usual proof of local existence of classical solutions, essentially going back to [1]. In this contribution we follow the presentation in [17]. With regard to the issue of higher regularity we have to use an additional induction loop.
Motivated by a preprint version of this paper the existence of global solutions for small initial data has already been proved in [6].
2. Heuristic derivations of models including radiation damping 2.1. The Vlasov-Poisson system with radiation damping. How could such a model look like? In [5,Thm. 1.4] it is shown that for suitable solutions of the Vlasov-Maxwell system which are isolated from incoming radiation in the non-relativistic limit c → ∞, c the velocity of light, the total rate of radiated energy is given by 2 3c 3 |D(t)| 2 where D is the dipole moment of the non-relativistic charge distribution, i.e. of the dipole moment of the solution of the Vlasov-Poisson system. This result yields a mathematical formulation and rigorous proof of Larmor's formula, see [9, (14.22)], in the case of Vlasov matter.
Therefore the goal is to implement an additional term into the Vlasov-Poisson system causing this loss of energy. As already suggested in [11,12] we modify the Vlasov equation by incorporating a small additional force term Here f + and f − give the distribution density of two species of charged particles where the mass of the particles is set to unity and the charge of the particles is +1 and −1, respectively. These densities are functions of time t ∈ R and phase space variables (x, p) ∈ R 3 × R 3 . The resulting charge density is given by ̺, x, p) dp and ̺ = ̺ + − ̺ − , E is the electrostatic field generated by this charge density, and finally D is the dipole moment of the charge distribution, We will refer to this set of equations (1)-(3) as the Vlasov-Poisson system with radiation damping. The additional term in (1) is a generalization of the radiation reaction force term used in particle models; see [9, (16.8)]. We also note that for this system the quantity is decreasing, more precisely one obtains the subscript S referring to the name "Schott"-energy under which this energy can be found in the literature. This decreasing of energy is usually attributed to the effect of radiation damping.

2.2.
Problems with the Cauchy problem. In order to formulate the Cauchy problem of this system we of course have to give initial values for the particle densities f ± , but we also have to fix data for the dipole moment. Let us for a moment assume that we are equipped with a smooth solution of (1)- (3) such that the spatial support of f ± is compact on every time slice with constant time. We define the bare mass by x, p) dp dx .
Then (1) yields mass conservation ∂ t M = 0 as well as charge conservation ∂ t ̺ ± + ∇ · j ± = 0 for both species, where (5) j ± = pf ± dp are the current densities. Using equation (1) we compute the first three time derivatives of D and finḋ where D [2] is defined by Introducing the new variable y(t) :=D(t) and abbreviating ε : where we have to supply initial data f •,± and y 0 ∈ R 3 . While the physical setting only delivers initial data for the particle distribution the question arises how to specify y 0 in order to get a solution which is physically reasonable. We shall discuss now several approaches to deal with this problem.

2.3.
Methods of reduction. The Vlasov-Poisson system with radiation damping serves as an approximation of a relativistic system valid up to order c −3 or ε. On account of that it seems to be acceptable to modify this system by terms which are formally of higher order in ε. In [11] there are several ways of reduction discussed, each of it consisting of two basic steps in certain combinations: on the one hand a time derivative will be replaced modulo higher order terms via the Vlasov equation, e.g.D = D [2] + O(ε) and on the other hand a time derivative will be cancelled by a substitution in the characteristic equation X ± = P ± , P ± = ±(E + ε ... D). The so-called 2+1 reductionX ± = X ± ,P ± = P ± ∓ εD andD ≈ D [2] leads to the Vlasov equation (6) ∂ t f ± + (p ± εD [2] ) · f ± ∇ x ± E · ∇ p f ± = 0 which is thoroughly studied in [12]. Using a 3+0 reduction scheme, replacing ... D byḊ [2] , we find the Vlasov-equation (7) ∂ t f ± + p · ∇ x f ± ± E + εḊ [2] · ∇ p f ± = 0 .

A reduction via geometric perturbation arguments.
Here we want to give an alternative formal derivation of equation (7) using an analogy to single particle models. First note that in case of a single particle the problem of the physical reasonable initial value is discussed since the beginning of the last century, see e.g. [9, chapter 16]. In [13] it has been observed that in particle models this problem has the structure of a geometric singular perturbation problem, and the "physical" dynamics are obtained on a center-like manifold of the full dynamics. For sake of discussion let us assume for a moment that the first equation in (SGPP ε ) is an ODE instead of a PDE and to simplify notations in this subsection f denotes the couple (f + , f − ). The same convention will be used for f • = (f •,+ , f •,− ). Then theorems from singular geometric perturbation theory, see [10], support us with invariant manifolds. More precisely, the manifold M 0 = {y = h 0 (f ) = D [2] } being invariant under the dynamics of (SGPP ε=0 ) persists for small ε > 0: For sufficiently small ε > 0 there is a manifold M ε = {y = h ε (f )} being invariant under the dynamics of (SGPP ε ) which remains close to M 0 in a suitable sense and carries those trajectories which are bounded in y.
Unfortunately, this theory does not apply in our case because we are dealing with a phase-space of infinite dimension; thus the proof of the existence of an invariant manifold is hard. For sake of discussion we shall take the existence of a smooth invariant manifold for granted and assume that it is given by for all times the solution exits. In order to find an appropriate approximation of (SGPP ε ) let us assume that we are furnished with smooth maps h ε defining the manifolds M ε and that we are able to expand in ε, We plug (8) and (9) into the second line of (SGPP ε ). Assuming that we are allowed to exchange the order of differentiation with respect to t and expansion in ε we compute Equating powers of ε we obtain Note that D [2] is a function of E and ̺ ± . But because E itself is defined by means of the Poisson integral with source ̺ it is a function of f ± , too. Therefore D [2] = D [2] (f ). Neglecting higher order terms we also end with the Vlasov-equation (7), which shall be discussed in the following. i

The reduced Vlasov-Poisson system with radiation damping, main results
3.1. Reformulation using principal values. First we have to remove an obstacle. On the one hand . On the other hand in [11, p. 3579] it is argued that (7) is a borderline case for the usual proof of local existence of classical solutions: In order to obtain a local existence theorem the orders of differentiability of the different unknowns should fit together properly. If it is supposed that the potential u of E is ktimes differentiable, then f and ̺ are k − 2 times differentiable and thus by solving the Poisson equation ∆U = 4π̺ two orders of differentiability should be gained which is not the case in spaces of pointwise differentiable functions.
We shall get rid of both of these problems by using the Vlasov-equation (7) and integration by parts one more time, but know we are forced to use principal values. To this end we additionally introduce the i On a very formal level hε(f ) is given by ∞ j=0 ε j D [2],(j) , where, using the Vlasov equation, all derivative w.r.t to time can be replaced by derivatives w.r.t to phase-space. notation E ± = −∇ ̺ ± (y) dy |x−y| ; then some elementary calculations, using the identity ∇ x |y − x| −1 = −∇ y |y − x| −1 , lead to Using charge conservation ∂ t ̺ ± + ∇ · j ± = 0 we e.g. compute Here id denotes the identity mapping in R 3 and for a suitable function j the operator H is defined by Note that the kernel K(z) = −3 z⊗z |z| 2 + id is bounded on R 3 \ {0}, is homogeneous of degree zero and satisfies |z|=1 K(z) dσ(z) = 0. Thus, using the Calderón-Zygmund inequality, H is well defined for smooth functions with compact support and for 1 < p < ∞ can be extended to a bounded linear operator mapping L p (R 3 ) to L p (R 3 ), see [20].
We call the following set of equations the "reduced Vlasov-Poisson system with radiation damping" where ̺ ± , j ± and H(j ± ) are defined according to (2), (5) and (10). To be precise we shall define what we mean by a classical solution. We adopt the formulation of the usual Vlasov-Poisson system given in [17].
is a classical solution of the reduced Vlasov-Poisson system with radiation damping on the interval J ⊂ R if the following holds: (i) The functions f ± are continuously differentiable with respect to all its variables.
(ii) The induced spatial densities ̺ ± and j ± , the induced electric field E and the induced damping term D [3] exist on J × R 3 and J, respectively, and are continuously differentiable. (iii) For every compact subinterval I ⊂ J the electric field E is bounded on I × R 3 . (iv) The functions f ± , E and D [3] satisfy the equations in (rVPRD ε ).

Local solutions, uniqueness and higer regularity.
It is the main result of this paper to establish local existence and uniqueness of classical solutions of this system. Furthermore we want to establish higher regularity of solutions of (rVPRD ε ) depending on the smoothness of the initial data. To the best knowledge also in the case ε = 0, which is he usual Vlasov-Poisson system, this question is only investigated in [15] and not published elsewhere. Surprisingly in that paper only smoothness in spatial directions is proved whereas the methods also allow for proving smoothness for all space-time directions. At last we want to establish certain bounds holding uniformly in 0 ≤ ε ≤ 1. These bounds are needed in [2] in order to prove the asymptotic approximation of solutions of the Vlasov-Maxwell systems by densities which are built from solution of (rVPRD ε ). ii We shall take our initial data from the following class of smooth functions with compact support. Fix an integer k and constants R 0 and S 0 . With regard to the initial data we require that (11) f Here ||·|| W ∞,k (R 3 ×R 3 ) denotes the usual Sobolev-norm. For further notations also in the following theorem, see the beginning of Section 4. Now we can state the main result of this paper.
The constantsT and M (T ) are only depending on k, R 0 and S 0 . In particular they are independent of ε.
The proof is elaborated in Section 4.

Global solutions.
Next we turn to the question of existence of global in time classical solutions. Motivated by a preprint version of this paper global solvability for small initial data has already been proved in [6]. One crucial ingredient of every proof of existence of global solution with unrestricted initial data is an a-priori bound on the second velocity momenta p 2 f ± dp dx. iii This bound is usually obtained by some energy conservation equation. We introduce the notation and computeḊ [1] = D [2] + εM D [3] for a solution of (rVPRD ε ), where M is the time-independent bare mass, see (4). We recallḊ [2] = D [3] and define the energy of the system (rVPRD ε ) by This energy is decreasing, more precisely we have d dtẼ S = −ε|D [2] (t)| 2 .
ii For the approximation we also need some correction terms of 'Darwin-order' iii Very recently in [7] a proof of the existence of global classical solutions of the Vlasov-Poisson system has been published, which does not rely on an a-priori bound of the second velocity momentum. It remains to investigate wether this proof can be applied to the system at hand.
Obviously, this energy is indefinite. This situation is well known from the Vlasov-Poisson equation in the gravitational case, where the potential energy is negative. Here the question is, wether or not it is possible to estimate |D [1] · D [2] | against one of the immanently positive terms in a suitable way. Of course we can estimate 2|D [1] · D [2] | ≤ |D [1] | 2 η + η|D [2] | 2 for all η > 0; therefore, it would be sufficient to bound |D [1] | 2 . Unfortunately, using the usual bounds for velocity momenta, see e.g. [17, Lemma 1.8], we only get the estimate , which is just not sufficient for our purposes. We also can't make use of the smallness of ε; in order to swallow the |D [2] | 2 term we have to choose η = εM 2 and end with the factor 1 M in front of the |D [1] | 2 term which is just cancelled by the square of ||f ± || 1/2 1 in (13). It is because of this, that the existence of global classical solutions seems to be out of reach at the moment. If we take the boundedness of the second velocity momenta for granted there are no obvious obstacles to prove the existence of global solutions following the lines of the proof of Lions and Perthame in [16]. iv

Proof of Theorem 2
Theorem 2 will be proved in two steps. After collecting some well known estimates in Subsection 4.1 we define a modification of the usual iteration scheme in Subsection 4.2 and establish certain bounds on the iterates. In Theorem 12 in subsection 4.3, we prove that the iterates converge to a solution of (rVPRD ε ) and that this solution is unique. Finally, in Subsection 4.4 we shall prove higher regularity of solutions, see Theorem 14. Together Theorem 12 and Theorem 14 give Theorem 2.
Throughout this section we fix a positive integer k and some constants R 0 and S 0 . Generic constants which may be computed from k, S 0 and R 0 are denoted by C, C 1 , C 2 and so on. Constants C may change from line to line. For 1 ≤ p ≤ ∞ we denote the spatial-space and the phase-space L p norms by ||·|| p . If we consider a pair of functions then ||·|| p gives the sum of the norms of the pair, e.g. ||f ± || p = ||f + || p + ||f − || p . For t ∈ R we denote by f (t) the function and in the same sense we use ̺(t), j(t) and E(t). By C(R n ) and C k c (R n ) we denote the space of k times continuously differentiable functions on R n , the subscript c indicates compactly supported functions.
For differentiable functions f , we denote (the vector of) partial derivatives with respect to t, x, p or combinations of it by ∂ t f , ∂ x f , ∂ p f , ∂ x,p f and so on. Vectors of partial derivatives of order l will be denoted by ∂ l x,p f , ∂ l t,x,p f and so on. The subscripts indicate which partial derivatives are included.
For ̺ ∈ C c (R 3 ) and j ∈ C c (R 3 ; R 3 ) we define D [3] ̺j : where c CZ is the Calderón-Zygmund constant. Let J ⊂ R an intervall, ̺ ∈ C(J; C c (R 3 )) as well as j ∈ C(J; C c (×R 3 ; R 3 )). Then the function J ∋ t → D [3] ̺j (t) := D [3] ̺(t)j(t) is continuous. If additionally ̺ and j are continuously differentiable with respect to t, then D [3] ̺j is differentiable and the usual product rule holds true.
For an intervall J ⊂ R and a mapping F ∈ C(J × R 3 ; Furthermore, we set For an initial datum f • according to (11)  Lemma 5. Assume that F is additionally continuously differentiable with respect to x and bounded on I × R 3 for every compact interval I ⊂ J. Then and for all (s, t) ∈ J × J the map Z F (s; t, ·, ·) is a measure preserving C 1 -diffeomorphism with inverse Z F (s; t, ·, ·) −1 = Z F (t; s, ·, ·).

4.2.
The iteration scheme and estimates on the iterates. We set up a modification of the usual iteration scheme. Fix initial data f •,± according to (11). For n = 0, s, t ∈ [0, ∞) and x, p ∈ R 3 we define Z ± 0 (s; t, x, p) = (x, p). If Z ± n is defined we set f ± n (t, x, p) = f •,± (Z ± n (0; t, x, p)) P n (t) = max P Z + n (t), P Z − n (t) X n (t) = R 0 + t 0 P n (s) ds ̺ ± n (t, x) = f ± n (t, x, p) dp j ± n (t, x) = pf ± n (t, x, p) dp Lemma 6. The iteration scheme is well defined and for all t ≥ 0, n ∈ N and x, p ∈ R 3 and ε ∈ [ 0 , 1 ] we have Proof. Using Proposition 3, Lemma 4 and Lemma 5 this Lemma follows immediately by induction on n ∈ N.
Next we shall give an estimate of the size of the support uniformly in n. Let P : [0, a) → R denote the maximal solution of the equation and define Lemma 7. For all 0 ≤ t < a, n ∈ N and ε ∈ [0, 1] we have P n (t) ≤ P (t) and X n (t) ≤ X(t) .
Proof. Using Lemma 5 (ii) and Lemma 6 (iv) and (v) this follows by induction on n ∈ N.
In the next steps we prove that a bound on the p-support on an interval J ⊂ [0, ∞), 0 ∈ J, implies bounds on all other quantities as well as the existence of a smooth solution of (rVPRD ε ) on the interval J. We assume that we are furnished with such an interval J and a continuous monotonously increasing positive function, called P again, such that for all n ∈ N, ε ∈ [0, 1], t ∈ J and x, p ∈ R 3 P n (t) ≤ P (t) .
In fact, this is at least true for J = [ 0 , a ). In addition we introduce the following helpful conventions. Henceforward, continuous monotonously increasing positive functions functions, which can be calculated by the knowledge of P and the basic constants k, S 0 and R 0 alone, are denoted by C(·), C 1 (·), C 2 (·). Especially they are independent of n ∈ N and ε ∈ [0, 1]. While C 1 (·), C 2 (·) will be fixed functions, the definition of C(·) may change from line to line, e.g. we have C(t) = P (t) = P (t) exp(C(t)) + C(t).
Proof. Using Lemma 7 and Lemma 5(iii) we compute Plugging this estimate together with Lemma 6(ii) in Proposition 3(iii) we have where we choose C 1 (·) in such a way that additionally We define C 2 (·) as the solution of By induction on n we conclude that ̺ ± n (t) ∞ ≤ C(t) for all t ∈ J and all n ∈ N.
Utilizing Lemma 8 the electric fields can be estimated by

Employing Lemma 4 we have for the dipols
Collecting these two estimates we conclude and using Gronwall's Lemma Combining (15) and (16) we have By induction it follows that Plugging this inequality into (16) we get For a compact intervall I ⊂ J we define Corollary 10. For every compact interval I ⊂ J the sequences (f ± n ) n , (Z ± n ) n , (̺ ± n ) n , (j ± n ) n , (E n ) n and (D [3] n ) n are Cauchy sequences uniformly in I × R 3 × R 3 , resp. ∆ I × R 3 × 3 , resp. I × R 3 , resp. I and ε ∈ [0, 1].
Proof. This is an immediate consequence of Lemma 9 in combination with Lemma 6.
Lemma 11. The sequence (∂ x E n ) n is Cauchy sequence uniformly in I × R 3 and ε ∈ [0, 1] for every compact interval I ⊂ J.
Proof. Using Proposition 3(iii) for every 0 < d ≤ R we have Choosing d sufficiently small proves the claim. 4.3. C 1 -solution and estimates. We define Theorem 12.
Remark 13. The usual continuation property also holds for (rVPRD ε ): an a-priori bound on the velocity support yields global existence of solution. This can be shown in the same way as in [17, Step 7, p. 401].
The crucial point is that the constants C 1 and C 2 in (14) do only depend on the L 1 and the L ∞ norm and not on the size of the support of the initial data.
4.4. C k -solutions and C k -estimates.
Theorem 14. f ± , ̺ ± , j ± , E and D [3] are C k and for all integers 0 ≤ j ≤ k we have the estimate [3] (t)| ≤ C(t). We shall proof this theorem by induction on 1 ≤ l ≤ k. Before starting we give a precise statement of the induction hypothesis.
Definition 15. For 1 ≤ l ≤ k we say H(l) holds true, iff the following holds: (i) For all 1 ≤ j ≤ l we have x,p f ± n ) n , (∂ j x,p Z ± n ) n and (∂ j x ̺ ± n ) n converge to ∂ j x,p f ± , ∂ j x,p Z ± and ∂ j x ̺ ± uniformly in ε ∈ [ 0 , 1 ] and uniformly on I × R 3 × R 3 , ∆ I × R 3 × R 3 and I × R 3 , respectively, for all compact subintervals I ⊂ J.
For the third step we apply Lemma 3(iii) on ∂ l x E n and ∂ l x ̺ n and chose d sufficiently small again. Thus, H(l + 1)(iii) is shown.