A local sensitivity and regularity analysis for the Vlasov-Poisson-Fokker-Planck system with multi-dimensional uncertainty and the spectral convergence of the stochastic Galerkin method

We study the Vlasov-Poisson-Fokker-Planck (VPFP) system with uncertainty and multiple scales. Here the uncertainty, modeled by multi-dimensional random variables, enters the system through the initial data, while the multiple scales lead the system to its high-field or parabolic regimes. We obtain a sharp decay rate of the solution to the global Maxwellian, which reveals that the VPFP system is decreasingly sensitive to the initial perturbation as the Knudsen number goes to zero. The sharp regularity estimates in terms of the Knudsen number lead to the stability of the generalized Polynomial Chaos stochastic Galerkin (gPC-SG) method. Based on the smoothness of the solution in the random space and the stability of the numerical method, we conclude the gPC-SG method has spectral accuracy uniform in the Knudsen number.

1. Introduction. In this paper, we are interested in the Vlasov-Poisson-Fokker-Planck (VPFP) system with multi-dimensional random inputs. The VPFP system describes the Brownian motion of a large system of particles in a surrounding bath, with a wide range of applications in plasma physics [4]. The physical system usually contains uncertainty, which can not be precisely described by deterministic partial differential equations. In this paper, we will mainly focus on the system with random initial input due to measurement errors or random impurity of the environment, and study how the randomness will affect the physical system. The uncertainty is modeled by multi-dimensional independent random variables with given probability density functions. We first study how the random initial data propagate in time, as well as the long-time behavior of the solution. We also study the stability and the convergence rate of the numerical method to the VPFP system with uncertainty, specifically, the generalized Polynomial Chaos stochastic Galerkin (gPC-SG) method. Both problems need an understanding of the regularity of the solution in the random space.
There are plenty of developments regarding the solution of elliptic or parabolic equations with uncertainty [2,5,6], while the regularity of the solution in the random space to kinetic equations has seldom studied until recently [9,15,13,16,19,17]. Kinetic equations give a uniform description of both mesoscopic and macroscopic physical quantities in terms of the Knudsen number . When → 0, the microscopic kinetic model approaches to the macroscopic hydrodynamic models [3]. Numerical and analytical difficulties increase when is small. In the numerical aspect, one efficient multiscale paradigm is the Asymptotic-Preserving schemes, which mimic the asymptotic transitions from kinetic equations to their diffusion or hydrodynamic limits in the numerically discrete space [12]. [14] extended this concept to random kinetic equations by stochastic Asymptotic-Preserving (sAP) methods. [13] gives the first result of uniform-in-regularity in the random space for linear transport equation, [15] gives the first result on uniform regularity for a nonlinear system, the Vlasov-Poisson-Fokker-Planck system. There are also some other uniform regularity results, we mention some of them here [16,22,10,17,19,17].
Depending on different scales, the VPFP system possesses two distinguished asymptotic limits, the high field limit, and the parabolic limit [1]. For the deterministic VPFP system without scaling parameters, [8,18,20] studied the convergence of the weak solution to its asymptotic limits, while [11] gave regularity results for classical solution near the global Maxwellian. For the VPFP system with uncertainty and scaling parameters, [15] get an exponential decay of the perturbative solution independent of the small parameter under some mild initial condition, which leads to a uniform regularity of the solution in the random space for both high field and parabolic limits.
This paper studies the same VPFP system as [15], but in more physically interesting setup. Space and velocity variables are in R 3 , and the random variable is d-dimensional. As to the long-time behavior, the first improvement compared to [15] is a much sharper decay rate to the global Maxwellian under a more general condition on initial data. The upper bound for the decay of the random perturbative solution near the global Maxwellian in a suitable Sobolev norm is O e −t/ in [15], which is improved to O e −t/ 1+a , where a = 1 for parabolic regime and a = 0 for high field regime, under the condition that the initial perturbation is smaller than O −(1+a) . This implies that for both regimes, as goes to 0, The range of the initial perturbation that will decay exponentially becomes larger. While in [15], only small random perturbation will decay exponentially in high field regime. Other than obtaining a sharp estimate when is small, getting rid of the bad dependency on large M in the initial condition for solutions in Sobolev norm H M z (L 2 x,v ) is another improvement in this paper. [15] mentioned briefly about eliminating the dependency on M for the one-dimensional case at Appendix, we generate it in detail for multi-dimension. Why a sharp estimation on H M z (L 2 x,v ) in terms of small and large M is not trivial and how we overcome the difficulties is stated in Section 2.3.
We treat the high field regime and the parabolic regime in a unified framework in this paper. The regularity of the solution in the random space comes from the study of the sensitivity of the perturbation near the global Maxwellian. With carefully designed energy norms, we combine the microscopic energy estimate and the macroscopic one to get the proper Lyapunov-type inequalities, which allows us to obtain the uniform regularity in the random space in terms of the scaling with initial data Here the distribution function f (t, x, v, z) depends on time t, position x, velocity v and random variable z. The density function ρ(t, x, z) is defined as, The collision operator F describes the Brownian motion of the particles, where M is the global Maxwellian, In the dimensionless system, δ represents the reciprocal of the scaled thermal velocity, while is the scaled thermal mean free path [20]. We will introduce two different regimes for this system [1]. One is the high field regime, where δ = 1.
Another is the parabolic regime, where δ = . We will study the two different scalings in a unified framework, where we assume, The random perturbation introduced by the initial data is characterized by a d-dimensional random variable z ∈ I z ⊂ R d . It is in a properly defined probability space (Σ, A, P), whose event space is Σ and is equipped with σ−algebra A and probability measure P. Define π(z) : I z → R + as the probability density function of the random variable z(ω), ω ∈ Σ. So one has a corresponding weighted L 2 space in the measure of dµ(z) = π(z)dz.
We further define according to the dependent variables of f and here · 2 2 is the regular Euclidean norm for vector function f ; and f , g µ is the corresponding inner product.

The perturbative solution.
It is easy to check that the global Maxwellian is a stationary solution to the VPFP system. We further introduce the perturbative solution h(t, x, v, z), perturbative density σ(t, x, v), perturbative flux u(t, x, z), near the the global Maxwellian, The perturbative solution h satisfies, where L is the so-called linearized Fokker-Planck operator, It is straightforward to see that multiplying 1, v to (4), then integrating it over v gives We call (4)-(5) the microscopic system, and (7)-(8) the macroscopic system. Moreover we define the projection onto the null space of L, N (L) = Span{ √ M} as, and one can check that, the random perturbation will exponentially decay in time. This assumption is not reasonable because stronger assumption is required on smoother initial random perturbation. Therefore, the estimate on h 2 H M z is not sharp enough for large M . We will overcome this difficulty by adding a weight to ∂ M z h as follows, where q β is defined as, here d is the dimension of random variable z. The weight q β is used to eliminate To sum up, we obtain a sharp estimate in both small and large M on then integrate the final result over dµ(z) to get Theorem 2.1. Before we present the main theorem on the sensitivity analysis, we first introduce some constants that will be used frequently later. Notation.
-C S : For ∀h ∈ H 2 x (L 2 v ), by Sobolev embedding for ∀h ∈ H r z (H 2 x ) with r the smallest constant strictly larger than d 2 , by Sobolev embedding, where q i is defined in (19). See Appendix A for the boundedness of the constant A. -C N : then the analytic solution (h, ∇ x φ) decays in time as follows, Here all the constants are independent of and M , where (10), (20), (22), (23) respectively. Remark 1. This theorem implies two things about the VPFP system with uncertainty.
1. The random perturbation around the steady state will exponentially decay. As → 0, the VPFP system becomes less sensitive to the random perturbation. 2. The regularity of the solution to the VPFP system in the random space is preserved. Furthermore, the regularity result is uniform in .
3. Sensitivity analysis for the solution with initial random perturbation.
In this section, we will prove Theorem 2.1. Theorem 2.1 is about the energy in , however, we will first prove the energy in H N x (L 2 v ) for {∂ β z h} |β|≤M , and then one just need to do another integral over For fixed z, we define the weighted energy as, where q β is defined in (19). Accordingly, we define the dissipation terms as, The main strategy is to do energy estimates on the following micro-macro systems. By taking ∂ β z , |β| ≤ M to (4) -(5) and (7) -(8), one has the micro-macro system, and, Combine the two energy estimates in Lemma 3.4 and 3.5 in a proper way, one can obtain a Lyapunov-type inequality as in Lemma 3.6 for which is exactly the energy we want to estimate. Finally apply the continuity argument to Lemma 3.6, one can obtain the sensitivity result in Theorem 2.1.
To get the optimal estimates, one needs to carefully deal with the two nonlinear terms in (26) and (27), that is where q, C S , A, C N are constants defined in (19), (20), (22), (23).
Proof. First note, where the first inequality comes from the Sobolev embedding (20) and While the last inequality is true for ∀N ≥ 3. Therefore, for e = 1−a 2 . One obtains the last inequality by using Young's inequality and = c 2 β , then the last second term of (33) can be bounded by |β|≤M i≤β One can do similar estimates to the last term. Therefore, applying Lemma 3.1 to (33) and breaking where q, C S , A, C N are constants defined in (19), (20), (22), (23).
Proof. We first sum over |α| ≤ N , Noticing by changing j and α − j, i and β − i, the following equality holds, Therefore, (36) shows that the first term and second term in (35) can be cancelled out when summing over |β| ≤ M , i ≤ β. Therefore one just needs to prove the bounds for We first sum over |α|≤N j≤α , similar to (31), one separates {j ≤ α} into |j| ≤ |α| 2 and |j| > |α| 2 , then using the Sobolev Embedding on ∂ j respectively. Then one ends up with the following bound, Afterwards, for |β|≤M i≤β , similar to (33)-(34), it is straightforward to get, The following is some equalities and inequalities that will be frequently used later.
Proof. If one takes ∂ α x and multiplies (26), then integrates it over x, v, sums over |α| ≤ N , |β| ≤ M , by Proposition 2 (a), Proposition 1 and Lemma 3.2, it is straightforward to get the above Lemma.
Proof. If one takes ∂ α x and multiplies (28), then integrates it over x, finally sums it over |α| ≤ N , |β| ≤ M , by Lemma 3.3, and Proposition 2 (c), (d), it is straightforward to get, Then by Proposition 2 (b), (e), one has, Let γ = 2 a , by the fact that a ≤ 1 a , one completes the proof.

The gPC-SG Method.
4.1. The numerical method. In this section, we will review a numerical method gPC-SG and apply to the VPFP system with uncertainty. We will prove the stability and the spectral accuracy of the method. For random variable z = (z 1 , · · · , z d ), if z i , 1 ≤ i ≤ d are independent of each other and the probability density function of z i is π i (z i ), then Therefore, let {Φ i k } ∞ k=0 be the corresponding orthogonal polynomial basis with respect to π i (z i )dz i , then the orthogonal polynomial basis for µ(z) can be written as, {Φ i } |i|≥0 satisfies the orthogonal condition under the measure µ(z), The K-th order subspace is therefore spanned by {Φ i } |i|≤K . As a popular numerical method, the generalized Polynomial Chaos stochastic Galerkin (gPC-SG) method is to find the approximate solution in the truncated K-th order subspace. Insert the approximate solution f K ,φ K in the form of, into (1) and then project it onto the subspace, one has the system for the approximate solution, where Φ K (z) = {Φ i } |i|≤K is the vector function that contains all basis functions up to K-th order. Similarly the approximation for the perturbative solution h(t, x, v, z) is defined as,ĥ and correspondingly the approximation for the perturbative density and flux, or equivalently, the deterministic coefficients ofĥ K , i.e. the vector function ĥK , |β|≤K satisfies the following deterministic system, where Theorem 4.1 and Theorem 4.2 about the approximate solution by gPC-SG are based on the following condition of the basis functions Φ β (z).
Remark 3. This is a generalization of the condition given in [19]. The i.i.d normalized Legendre polynomials, which corresponds to Uniform distribution in [−1, 1] d with pdf π(z i ) = 1 2 d , and the Chebyshev polynomials, which corresponds to the random variable with pdf π(z i ) =

4.2.
Main results on stability and accuracy of the gPC-SG method. We want to estimate the error of the approximationĥ K obtained by the gPC-SG method. We first decomposed the error into two parts, The first part of the error h K p is caused by the gPC projection, which is related to the smoothness of the solution in the random space, which has been studied in Section 2 and 3. On the other hand, the second part of the error h K e is caused by the stochastic Galerkin, which is related to the stability of the gPC-SG method.
The difficulty in the proof of the stability of the gPC-SG method is to get a sharp estimate on ĥK 2 in terms of large K. If we directly estimate ĥK 2 , then similar to the sensitivity analysis in (15), the nonlinear term on the RHS of (51) will be as large as so only when the initial data is as small as the gPC method is stable. Actually, there is a much sharper estimates in terms of large K. Under Condition 1, we introduce a weight toĥ β as follows, Here the weight c β is used to eliminate max |κ|,|γ|,|β|≤K E κγβ −1 , while q β is used to eliminate K −1 in (54). Before we come to the final results on the stability and accuracy of the gPC-SG method, we will first define some constants that will be frequently used later. Notation. -Â:Â with C 2 defined in (53), q defined in (19).
Theorem 4.1. (Stability of the gPC-SG method) Let under Condition 1, for ∀N ≥ 3, if initially, then the approximation solution (ĥ K , ∇ xφ K ) obtained by gPC-SG method decays in time as follows, Here all the constants are independent of and K, whereÂ, ξ, C N , C S are the same constants in Thoerem 2.1.

Remark 4.
The above theorem tells us that the gPC-SG method is stable under a mild condition on the initial randomness.

Remark 5.
Notice there is another sufficient initial condition directly on (h(0, x, v, z), ∇ x φ(0, x, z)) to guarantee the stability of the gPC-SG method, which has been derived in Remark 7.
Based on the regularity of the solution in the random space as in Theorem 2.1 and the stability of the gPC-SG method as in Theorem 4.1, one has the following Theorem about the spectral convergence of the gPC-SG method.
, then the error decays in time as follows, Here r is the smallest integer strictly larger than d 2 , Remark 6. The above Theorem tells us as long as the initial data (h, ∇ x φ) has enough regularity in the random space, and the initial perturbation around the global Maxwellian is smaller than O( 1 1+a ), then the gPC-SG method enjoys spectral accuracy.

5.
Stability of the gPC-SG method. In this section, we will prove Theorem 4.1. We will study the stability of the gPC-SG method in terms ofÊ K,N h ,Ê K,N ∇xφ defined as follows, whereĥ β is the β-th component of the vector functionĥ K (t, x, v), same forφ β . Alsô D K,N h1 ,D K,N σ ,D K,N ∇xφ are the corresponding dissipations of (1 − Π)ĥ K ,σ K , ∇ xφ K in the norm of · H N ω . One needs to do energy estimation on (51) to get the estimates forÊ K,N h andÊ K,N ∇xφ . Actually, if one compares (51) with (26), one finds the only difference is the nonlinear term. Therefore, one can use the same proof strategy to , except we need to re-estimate the nonlinear term again. Before we estimate the nonlinear term, we first define a characterized function χ γκβ , which implies, Here the first inequality comes from the orthogonality of the basis. The second one is because of the first property of Φ γ L ∞ z in Condition 1. Moreover, from the definition of χ, one has the following Lemma.
Compare Lemma 5.2 with Lemma 3.2 and 3.3, one notes that the estimates for the two nonlinear terms are similar, so one ends up with the similar energy estimates forÊ K,N in Theorem 4.1.

Remark 7.
Here we derive a sufficient condition on the initial data (h(0, x, v, z), ∇ x φ(0, x, z)). We requireÊ in whereĥ β (0) = h(0)Φ β dµ(z). By the classical approximation theorem, we know that for h ∈ H M,N , for some constant D depending on the measure µ(z). Plug (67) into (68), one gets, Similarly, one can get the bound for ∇ xφβ , Therefore the condition (66) becomes, hence a sufficient initial condition for the stability of the gPC-SG method is, If c β grows algebraically, then there exists an positive integer p, such that, c β (|β| + 1) p up to a constant, which implies that as long as M > d 2 + q + p, then To sum up, another sufficient condition to enjoy the stability of gPC-SG method as in (25) is that firstly the bound of the basis Φ β grows algebraically as Φ β L ∞ z (|β| + 1) p , and initially for N ≥ 3, M > d 2 + q + p.
6. Spectral convergence of the gPC-SG method. In this section, we will prove Theorem 4.2. Let us define the projection of the analytic perturbative solution h onto the subspace Φ K = {Φ β } |β|≤K as, Then we can decompose the error of the approximation solution ĥK , ∇φ K in the subspace Φ K into two parts, where h K p , ∇ x φ K p represents for the projection error, h K e , ∇ x φ K e are errors from the stochastic Gelarkin. Define vector functions, energy and dissipation terms for error as follows, where ∇ x φ e,β = ∇ x φ −φ K Φ β dµ(z); Here is the proof sketch of the spectral convergence of the gPC-SG method. Project the microscopic system for the perturbative solution (4)-(5) onto the truncated subspace {Φ K }, and then subtract the approximate perturbative system (50) from it, one has the following microscopic error system, If one does energy estimate to the above system, one has microscopic error estimates as in Lemma 6.1. If one multiplies v √ M, and integrate over v to (72), then one has the following macroscopic system of error, If one does energy estimate to it, one will obtain estimates as Lemma 6.2. If one combines the microscopic and macroscopic error estimates in a proper way as in Lemma 6.3, and then based on Corollary 1, one can obtain the spectral convergence of the gPC-SG method. From Theorems 2.1 and 4.1, one can derive the following Corollary.
Proof. If one multiplies h e,β to the β−th row of (72), integrates it over x, v, and sums over |β| ≤ K, one has, By (30) and integral by parts, the nonlinear term becomes, RHS of (79) = − 1 First notice that For I, one can use the definition of ω β in (55) and the bound for E γκβ in (62), bound it in a similar way as in Lemma 5.2.