Stability of a point charge for the Vlasov-Poisson system: the radial case

We consider the Vlasov-Poisson system with initial data a small, radial, absolutely continuous perturbation of a point charge. We show that the solution is global and disperses to infinity via a modified scattering along trajectories of the linearized flow. This is done by an exact integration of the linearized equation, followed by the analysis of the perturbed Hamiltonian equation in action-angle coordinates.


Vlasov-Poisson near a point charge
This article is devoted to the study of the time evolution and asymptotic behavior of a three dimensional gas of charged particles (a plasma) that interact with a point charge. Under suitable assumptions this system can be described via a measure M on R 3 x × R 3 v that is transported by the long-range electrostatic force field created by the charge distribution, resulting in the Vlasov-Poisson system Since this equation is rotationally invariant, the Dirac mass M eq = δ = δ (0,0) (x, v) is a formal stationary solution and we propose to investigate its stability. We consider initial data 1 of the form M = q c δ + q g µ 2 0 dxdv, where q c > 0 is the charge of the Dirac mass and q g > 0 is the charge per particle of the gas, which results in purely repulsive interactions. We track the singular and the absolutely continuous parts of a solution as M (t) = q c δ (x(t),v(t)) + q g µ 2 (t)dxdv, which yields the coupled system where λ = q 2 g /( 0 m g ) > 0, q = q c q g /(2π 0 m g ) > 0, q = q c q g /( 0 m c ) > 0 are positive constants 2 . 1.1. Main result. Our main result concerns (1.2) with radial initial data, where the point charge is located at the origin. For sufficiently small initial charge distributions µ we establish the existence and uniqueness of global, strong solutions and we describe their asymptotic behavior as a modified scattering dynamic. While our full result can be most adequately stated in more adapted "action-angle" variables (see Theorem 1.6 below on page 6), for the sake of readability we begin here by giving a (weaker) version in standard Cartesian coordinates: Theorem 1.1. Given any radial initial data µ 0 ∈ C 1 c (R * + × R), there exists ε * > 0 such that for any 0 < ε < ε * , there exists a unique global strong solution of (1.2) with initial data (x(t = 0),v(t = 0)) = (0, 0), µ(t = 0) = εµ 0 . 1 Here the initial continuous density f 0 = µ 2 0 is assumed to be non-negative, a condition which is then propagated by the flow and allows us to work with functions µ in an L 2 framework rather than a general non-negative function f in L 1 -see also our previous work [18] for more on this. 2 In these formulas, 0 is the vacuum permittivity, mg the inertia of a gas particle and mc the inertia of the point charge. 1 Moreover, the electric field decays pointwise and there exists an asymptotic profile γ ∞ ∈ L 2 (R * + × R) and a Lagrangian map (R, V) such that µ(R, V, t) → γ ∞ (r, v), t → ∞.
(1) Our main theorem is in fact much more precise and requires fewer assumptions, but is better stated in adapted "action angle" variables. We refer to Theorem 1.6.
(2) The Lagrangian map can be written in terms of an asymptotic "electric field profile" E ∞ : The first term corresponds to conservation of the energy along trajectories, the second term comes from a linear correction and the third term on the first line comes from a nonlinear correction to the position. This can be compared with the asymptotic behavior close to vacuum in [18,28] by setting q = 0.
1.1.1. Prior work. In the absence of a point charge, the Vlasov-Poisson system has been extensively studied and we only refer to [1,22,30,32] for early references on global wellposedness and dispersion analysis, to [7,18,28] for more recent results describing the asymptotic behavior, to [12,31] for book references and to [2] for a historical review. The presence of a point charge introduces singular electric fields and significantly complicates the analysis. Nevertheless, global existence and uniqueness of strong solutions when the support of the density is separated from the point charge has been established in [24], see also [4] and references therein, while global existence of weak solutions for more general support was proved in [9] with subsequent improvements in [20,21,25]. We also refer to [8] where "Lagrangian solutions" are studied and to [5,6] for works in the case of attractive interactions. Concentration, creation of a point charge and subsequent lack of uniqueness were studied in a related system for ions in 1d, see [23,33]. To the best of our knowledge, there are no works concerning the asymptotic behavior of such solutions.
The stability of other equilibriums has been considered for the Vlasov-Poisson system with repulsive interactions, most notably in connection to Landau damping [3,11,15,27]. In the case of Vlasov-Poisson with attractive interactions, there are many more equilibriums and their linear and nonlinear (in)stability have been studied [13,14,19,26,29], but the analysis of asymptotic stability is very challenging. We also refer to [17] which studies the stability of a Dirac mass in the context of the 2d Euler equation.
1.1.2. Our approach. In previous works on (1.2), the Lagrangian approach allows to integrate the solutions against characteristics but faces the problem of a singular electric field, while a purely Eulerian method leads to a poor control of the solutions, which makes it difficult to study the asymptotic behavior. In this paper, we introduce a different method based on the decomposition of the Hamiltonian to rewrite (1.2) as where the linearized Hamiltonian H 0 is given in (2.1) and the nonlinear Hamiltonian H pert corresponds to the self-generated electrostatic potential (1.17). In short, our approach combines a Lagrangian analysis of the linearized problem with an Eulerian PDE framework in the nonlinear analysis, all the while respecting the symplectic structure. This amounts to considering solutions as superpositions of measures on each trajectory of the linearized flow instead of measures on the whole phase space.
On a technical level, one faces the two difficulties of a singular transport field and the nonlinearity separately: the singular electric field created by the point charge is present in the linearized equation coming from H 0 , which is integrated exactly. The nonlinearity comes from the perturbed Hamiltonian H pert , but this leads to a simple nonlinear equation, with a nonlinearity which is smoothing.
More precisely, in a first step we analyze the characteristic equations of the linear problem associated to (1.2). These turn out to be the classical ODEs of the Kepler problem, which can be integrated in adapted "action-angle" coordinates. In these, the geometry of the characteristic curves is straightened and the linear flow is solved explicitly as a linear map. To treat the nonlinear problem, we conjugate by the linear flow and study the resulting unknown in an Eulerian, L 2 based PDE framework, based on energy estimates as in our recent work on the vacuum case [18]. This allows us to propagate the required regularity and moments to obtain a global strong solution. Moreover, we can readily identify the asymptotic dynamic: in a mixing type mechanism, the dependence on the "angles" is eliminated from the asymptotic electrostatic fields, and the scattering of solutions is modified by a field defined in terms of the "actions".
We remark on some features and context of our techniques.
(1) Since the system (1.2) is Hamiltonian and we solve the linearized system through a canonical change of unknown (i.e. a diffeomorphism respecting the symplectic structure), the nonlinear problem becomes quite simple after conjugation, see (1.16). (2) The moments we propagate are conserved by the linearized flow, unlike the physical moments in r , v . In fact, even in the nonlinear problem it is quite direct to globally propagate moments in action-angle variables, which already gives the existence of global weak solutions. (3) The asymptotic dynamic is easy to exhibit in action-angle variables through inspection of the formulas for the asymptotic electrostatic fields (see (1.18)). (4) It is notable that we do not require any separation between the point charge and the continuous distribution µ, addressing a question raised in [9, p. 376], (see also [24]). (5) We expect the methods presented here, based on integration of the linearized equation through "action-angle" coordinates, to be broadly applicable, both for local existence of rough solutions and especially for the analysis of long time behavior whenever the linearized equation corresponds to a completely integrable ODE without closed trajectories. This should include a large number of radial problems for plasmas since 1+1 Hamiltonian ODEs can be integrated by phase portrait. (6) The usefulness of action-angle variables for the Vlasov-Poisson equation was already exhibited in [10,13,16] where the authors produce a large class of 1d BGK-type waves which are linearly stable.
1.1.3. Remarks on the physical setup. Our primary interest here is to study the interaction of a gas of ions or electrons interacting with a (similarly) charged particle, subject only to electrostatic forces. In this case, up to rescaling, we may assume that λ = 1 in (1.2). Taking into account gravitational effects, we may also consider the more general case of a large charged and massive point particle with mass m c and charge q c interacting with a gas of small particles with mass-per-particle m g and charge-per-particle q g subject to both gravitational and electrostatic forces. In this case, our result holds whenever the principal gas-point particle interaction is repulsive, i.e. when (in appropriates physical units) whereas (due to the small data assumption on the gas at initial time) the gas-gas interaction may be repulsive (λ = 1) or attractive (λ = −1) depending on the sign of (q g ) 2 − (m g ) 2 .
The situation that is beyond the scope of our analysis is the case when the inequality in (1.3) is reversed and some trajectories of the linearized system are closed. Note that in this case, even the local existence theory is incomplete.

1.2.
Overview and ideas of the proof. We note that in the particular case of radial initial conditions 3 µ(x, v, t = 0) = µ 0 (|x| , |v|), (x,v)(t = 0) = (0, 0), (1.4) the Dirac mass in (1.2) does not move (i.e.x(t) =v(t) = 0) and the continuous particle distribution µ(t) of the solution is a radial function. The equations (1.2) reduce to the following system for µ(x, v, t): Per a slight abuse of notation with r := |x|, (r, t) = (x, t), by radiality the electric field E := −∇ x ψ of the ensemble can be computed as The "radial" phase space. Since as discussed the equations (1.2) are invariant under spherical symmetry and we will work with spherically symmetric data, it is more convenient to work on the . Note that r 2 v 2 drdv is the natural measure corresponding to that of radially symmetric functions on R 3 × R 3 , and hence we will work with the new density µ(r, v, t) := rvµ(r, v, t).
(1.7) This is chosen such that the (conserved) mass is the square of the L 2 norm of both µ on R * Moreover, the equations for µ simply read This equation is Hamiltonian and can be equivalently written as which leads to the conservation of energy The linearized system. In order to study (1.9), we first consider the linearized equation for a function f (r, v, t): This linear transport equation can be solved directly via its characteristic equationṡ One recognizes here the classical Kepler problem in the radial setting, which can be integrated using generalized 4 "action angle" coordinates (θ, a) (see also Figure 1). There exists a canonical transformation to "action-angle" coordinates: solves the free streaming equation (1.14) 1.2.3. Nonlinear analysis. Since (1.14) can be solved directly, we conjugate by this change of variables to stabilize the linearized system, thus defining γ as follows: γ(θ, a, t) := µ(R(θ + ta, a), V (θ + ta, a), t), Since the change of variable preserves the symplectic structure, we find that the full nonlinear problem (1.9) is equivalent to (1.16) Here the potential can be expressed in action angle coordinates as follows: Remark 1.4. While the trajectories of the linear equation (2.2) are straight lines in action-angle variables, in physical variables they correspond to an incoming ray followed by an outgoing one traced at varying velocities.
For the nonlinear problem this creates extra challenges, as interactions can occur over vastly disparate spatial scales. As the below Figure 2 illustrates, in some regimes the evolution R(θ + ta, a) is not a simple function of (θ, a), from which one of the two variables can be recovered once the other is known (see also Lemma 2.5 below).
It remains to study solutions to (1.16). The dispersion mechanism is accounted for through the conjugation with the linearized flow and we hope to show that this picture remains true when we add the remaining nonlinear contribution, i.e. we expect that solutions to (1.16) do not change too much  over time. We first use a bootstrap argument to propagate strong norms, which suffices to obtain global existence and decay of the electric field: Proposition 1.5. There exists ε * such that for all 0 < ε 0 ≤ ε 1 ≤ δ < ε * , the following holds. Let γ be a solution to (1.16) with initial data γ 0 on 0 ≤ t ≤ T and assume that for This in turn allows us to investigate the asymptotic behavior. It can easily be formally deduced once one observes that, given the bounds propagated by the bootstrap, on the support of the density, one has so that one expects that Ψ is asymptotically independent of θ: As a consequence, (1.16) becomes a perturbation of a shear equation: which can easily be integrated. This can be made rigorous under appropriate assumptions on the initial data and it leads to our main result.
Theorem 1.6. There exists ε 0 > 0 such that any initial data γ 0 satisfying We expect that the number of moments in (1.19) can be significantly reduced. It is interesting that the decay of the electric field can be obtained under much weaker assumptions (see Lemma 3.1), but it is unclear to us how much asymptotic information can be recovered in this case.
1.3. Organization of the paper. In Section 2, we study the ODE associated to the linearized equation and establish a number of geometric results and bounds on relevant quantities for the nonlinear problem. In Section 3, we study the nonlinear equation; we establish the moment bootstrap for weak solutions in Section 3.1, the derivative bootstrap for strong solutions in Section 3.1.2 and prove the modified scattering in Section 3.2 by first obtaining a weak-strong limit for scattering data in Section 3.2.1 and finally the convergence of the particle density in Section 3.2.2.
We emphasize that one can obtain decay of moments for weak solutions in a self-contained way using only the results of the Sections 2.1, 2.2 and 3.1.1.

Linearized equation
The goal of this section is to integrate the linearized problem (2.2) and to prove various estimates for the corresponding transfer functions. Lemma 1.3 follows easily from Lemma 2.1 below.

2.1.
Straightening the linear characteristics. The linearization of (1.10) at µ = 0 is the Hamiltonian differential equation associated to This is now a linear transport equation, which can be integrated easily once we know the trajectories of the corresponding ODE:ṙ where G : (1, ∞) → R and H : R + → (1, ∞) satisfy (2.5) These functions and related ones are studied in more details in Section 2.1.2. We can now solve the linear problem (2.2) via a canonical change of variable: in (2.4) defines a canonical diffeomorphism of the phase space (r, v) (with inverse (R(θ, a), V (θ, a)) as in (2.4)), which linearizes the flow in the sense that for the flow map Φ t (r, v) associated to the Hamiltonian Proof. That (Θ, A) and (R, V ) are inverse can be checked directly once one observes that r min is consistent: r min (A(r, v)) = r min (r, v) and r min (a) = r min (R(a, θ), V (a, θ)). It is direct to check that A is conserved along the flow. Moreover, ).
The first term on the right hand side gives A, and the second vanishes since when v = 0, G(r/r min ) = G(1) = 0, while direct computations show that the bracket in the last term vanishes. In addition, the same computations show that which shows that the transformation preserves the symplectic form and hence has Jacobian 1.

2.1.2.
Study of the structure functions. In this subsection, we study the geometric functions that arise from the change of variable in Lemma 2.1. These are independent of assumptions on the solutions.
Lemma 2.2. The functions G and H are almost linear In addition, we note the asymptotic behavior of G and its inverse Proof of Lemma 2.2. Since G can be integrated explicitly, we easily obtain the first line in (2.8) from (2.5), and the second line follows from the fact that H = G −1 . Now (2.9) follows by expanding the expression for G.
In addition, we will frequently consider first and second order derivatives.
Lemma 2.3. We have explicit formulas for the first order derivatives and the formulas for the second order derivatives , |∂ a ∂ a R(θ, a)| ≤ q a 4 ln a 2 q R(θ, a) .
Proof of Lemma 2.3. The formulas follow from direct calculations using that together with Lemma 2.2 to control the behavior of H.

Kinematics of linear trajectories.
Here we collect a few estimates on the behavior of the trajectories of (2.3). Using (2.7), we see that the linearized flow is simple in the action-angle variables. For simplicity of notation, given a function F (θ, a) in phase space, we will denote by F (θ, a) := F (θ + ta, a) (2.10) its evolution under the linear flow. This is a slight abuse of notation since the transformation depends on time; however all our estimates will be instantaneous so this should not lead to confusion. Since we will show that in action-angle variables, the new density is (almost) stable, we expect that the main role will be played by trajectories starting from the "bulk region" B defined by We start with a few simple observations. By definition, we have a universal lower bound for R, but in the bulk region, one can be more precise.
Lemma 2.4. We have the following control on R: (2.14) In addition, we have a more precise control in the bulk: when (θ, a) ∈ B, we have that |∂ a ∂ a R| a −4 ln ta 3 , and in particular, the change of variable a → R(θ, a) is well behaved.
Since the interaction involves quantities defined in the physical space, it will be useful to understand how to relate them to phase space variables. The next lemma is concerned with solutions of the equation for fixed t and r.

Lemma 2.5. Let
A := q r (1 + ), := c · min{1, r 3 for some fixed small constant c > 0 and define the regions In the region R 0 , we see that for any θ, there exists at most one a = ℵ(θ; r, t) solution of (2.15). In addition, we have that In the region R j , j ∈ {1, 2}, for each choice of a, there exists exactly one θ = τ j (a; r, t) solution of (2.15). In addition, we see that
We start with R 0 and denote a * = q/r, ϑ = θ + ta and x = a 2 |ϑ|/q so that we are considering the equation Now let a = a * √ 1 + h 2 for some h ≥ 0. We have that by (2.9) and Lemma 2.3 ds . (2.20) From this we see that if h is small enough, 0 ≤ h ≤ c, there exists a unique solution x to E 1 = 0, and this solution satisfies h ≤ x ≤ 3h. In addition, if h ≤ c/(q(a * ) 3 t), there holds that ∂ a R −1/(q(a * ) 3 ). Now in region R j , j ∈ {1, 2}, we see that and in particular, using (2.19), we find that In addition, we have that and the other bounds in (2.17) follow directly from the definitions. Finally, the last statement in (2.18) follows from the formula for ∂ a R in (2.20): Note that ϑ/ |ϑ| = 1 and letting x = b the bound where the term in parenthesis vanishes, then when 0 ≤ x ≤ 2b, there holds that q ≤ a 2 r ≤ H(2b), while for x ≥ 2b we see that both terms have same sign and H ≥ 0, so that ∂ a R(θ, a) ≥ tH (2b).

2.2.
Electric field and potential. Given an instantaneous density distribution µ(r, v, t), it is useful to introduce the "physical potential" Ψ of the associated electric field as in (1.9), explicitly given as when γ(ϑ, α, t) is the corresponding density distribution in action-angle coordinates as in (1.15). Then (1.17) can be rewritten as Ψ(θ, a, t) = Ψ( R(θ, a), t). This allows to obtain formulas for the derivatives of Ψ in action angle variables in terms of the electric field E, the local mass m and the density (compare also (1.9)): (r, t) := ∂ r m(r, t) = δ(R(ϑ + tα, α) − r) · γ 2 (ϑ, α, t)dϑdα.

(2.22)
Then for β ∈ {a, θ} we have (2.23) We note that the local mass has a trivial uniform bound but this can be made more precise.
Lemma 2.6. We can decompose m as m(r, t) = m s (r, t) + m n (r, t), where we have that for any , κ > 0 (2.25) Proof of Lemma 2.6. The decomposition corresponds to localizing in and out of the bulk zone defined in (2.11). Thus Using Lemma 2.4, we see that while using (2.12), We will use the following consequences: Proposition 2.7. There holds that (2.27) Proof of Proposition 2.7. Using (2.12) and (2.13), we find that .

2.2.1.
Study of the density. Controlling derivatives of γ requires estimates on the density; these are obtained in a similar way to the mass (see Lemma 2.6), but are more involved.
Lemma 2.8. The density can be decomposed into two terms, (2.28) The key observation is that the estimate for s only involves at most one copy of the large term a∂ a γ L 2 θ,a . Proof of Lemma 2.8. We can decompose into two regions m(r, t) = m 1 (r, t) + m 2 (r, t), where ϕ ≥1 (x) denotes a smooth function supported on {x ≥ 1/10} and equal to 1 for x ≥ 1/2. Study of ∂ r m 1 . This contains the main term. Integrating by parts, we observe that From Lemma 2.4 we recall that in the bulk region R ∼ at, and thus while on the other hand, using (2.12) and (2.11),
Remark 2.9. As can be seen from the proof of 2.8, we only need positive moments in a to control the area outside of the bulk where |θ| ≥ at. Such moments in a could be replaced by moments in θ, and thus positive weights in a are not necessary for our result.
Proposition 2.10. There holds that

30)
and (2.32) Proof of Proposition 2.10. The most important term is the term with mixed derivative (see Section 3.1.2). We recall from (2.23) that On the one hand, using Lemma 2.4, we see that and this leads to an acceptable contribution using Lemma 2.6 and Lemma 2.8. We now turn to Using Lemma 2.4 and (2.12), we see that and this term can be handled as before using Lemma 2.6 and Lemma 2.8. Finally, we compute that Using (2.13) and (2.14), we find that (1 + a 2 ) 1 + R(θ, a) + t 2 / R and that Using Lemma 2.6, we find that θ,a while using Lemma 2.8, we find that This finishes the proof.

Nonlinear analysis
In this section we consider the full nonlinear equation (1.16), We first establish global existence of strong solutions via a bootstrap in Section 3.1, then we demonstrate the modified scattering asymptotics in Section 3.2. This establishes Proposition 1.5 and Theorem 1.6.
3.1. Bootstraps and global existence. We first propagate global bounds using energy estimates.
The key property we will use is that the integral of a Poisson bracket vanishes. Commuting with appropriate operators gives the equations The key in the bootstrap estimates is that one can propagate a moments easily and that the terms with slowest decay involve only these moments (see m s in (2.25) and s in (2.28)). Interestingly, we will see in Section 3.1.1 that the moments can be bootstrapped on their own, allowing global bounds on weak solutions. These moment bounds allow to propagate another bootstrap for higher regularity. For simplicity, we only propagate the first order derivatives in Section 3.1.2.
3.1.1. Moment Bootstrap. It turns out that control of the moments can be bootstrapped independently of any derivative bound. Lemma 3.1. Let p ≥ 2 and assume that γ solves (1.16) on an interval 0 ≤ t ≤ T and assume that then there holds that Proof of Lemma 3.1. The moments can be readily estimated. By (3.2) we have that, for q ∈ R Using (2.26), the bootstrap assumptions (3.3) and Gronwall inequality, we find that Similarly, for q ≥ 0, using (3.2) and (2.27), we find that and we can again apply Gronwall estimate.
In particular, the case of ω = 1 gives the result of Proposition 1.5.
3.2. Asymptotic behavior. The analysis in this section is partially inspired by [18,28].
3.2.1. Weak-strong limit and asymptotic electric field. Before we obtain strong convergence of the particle distribution, we first need weak convergence of "asymptotic functions" which are defined in terms of averages along linearized trajectories. Given a bounded measurable function τ , we define τ (t) := τ (α)γ 2 (ϑ, α, t)dϑdα.
Integrating the time derivative then gives the bound (3.11).
The convergence of the scattering data allows to define the asymptotic electric potential and electric field