Particle simulation methods for the Landau-Fokker-Planck equation with uncertain data

The design of particle simulation methods for collisional plasma physics has always represented a challenge due to the unbounded total collisional cross section, which prevents a natural extension of the classical Direct Simulation Monte Carlo (DSMC) method devised for the Boltzmann equation. One way to overcome this problem is to consider the design of Monte Carlo algorithms that are robust in the so-called grazing collision limit. In the first part of this manuscript, we will focus on the construction of collision algorithms for the Landau-Fokker-Planck equation based on the grazing collision asymptotics and which avoids the use of iterative solvers. Subsequently, we discuss problems involving uncertainties and show how to develop a stochastic Galerkin projection of the particle dynamics which permits to recover spectral accuracy for smooth solutions in the random space. Several classical numerical tests are reported to validate the present approach.


Introduction
The Landau-Fokker-Planck kinetic equation, also referred to as the Landau equation, holds significant importance as a fundamental tool for understanding the complex dynamics of charged particles within a collisional plasma.The equation plays a crucial role in elucidating the behaviour of particles that interact over long distances, primarily influenced by Coulomb forces.These long-range interactions give rise to small-angle collisions between the particles, leading to a breakdown in the finiteness of the classical Boltzmann collision operator [4,15,26,40].
Due to its profound implications in the development of fusion reactors, the formulation of this model has had a far-reaching impact on both the theoretical and applied scientific communities.Researchers and engineers in the field have extensively utilized this formalized model to gain insights into the behaviour of charged particles in plasmas, paving the way for advancements in fusion reactor design and other related applications.
The construction of numerical methods for the Landau-Fokker-Planck equation has to deal with several challenges, among which are the high dimensionality of the problem, the structural properties, and the presence of multiple scales.As in most kinetic equations, the construction of numerical methods for plasma physics can be subdivided into two main categories, the first based on direct discretizations of the collisional PDEs, like spectral methods, see e.g.[13,18,19,23,33,34], and the second based on the approximation of the underlying particles' dynamics.In this latter direction we mention particle-in-cell methods [1], direct simulation Monte Carlo methods [2,3,5,6,35,39] and deterministic particle methods [7,8].It is worth to remark that Monte Carlo approaches based on the groundbreaking contribution of Bobylev and Nanbu [2,30] are consistent with the Landau equation thanks to their link with the Boltzmann equation when collisions become grazing.
To enhance prediction accuracy, it is crucial to understand how sensitive physical models are to potential uncertainties in the constitutive parameters that characterize the system behaviour and influences of the different scales.Consequently, addressing the issue of quantifying these uncertainties becomes a major concern, both analytically and computationally.The development of numerical techniques must consider the added complexity arising from the increased dimensionality of the system due to the presence of random inputs in the model.In this direction, several research efforts focused on the design of efficient numerical techniques for uncertainty quantification (UQ), see [16,20,21,25,31,37,38,44] and the references therein.We mention also recent results on UQ in plasma physics obtained in [12,22,24,42].
Among the most popular UQ techniques, stochastic Galerkin (sG) methods based on generalized polynomial chaos expansions have shown the ability to achieve spectral accuracy in the random space for smooth solutions [43].However, in contrast to stochastic sampling methods, these methods are highly intrusive and require to design new algorithms.Furthermore, a direct application of sG methods may lead to the loss of important structural properties, like conservation of physical quantities, nonnegativity of the solution, and hyperbolicity close to fluid regimes, see e.g.[14,17].To overcome this issue, a different approach based on projecting the particle dynamics through generalized polynomial chaos was pro-posed in [9,10,28,36].Recently, the methodology has been extended to plasma simulations with BGK collisions in [27].
In this article, our focus lies on the more realistic scenario where plasma collisions are described by the Landau equation, taking into account uncertain quantities.Specifically, we concentrate on the space homogeneous case, where the collision algorithm follows the approach outlined in [2,5] whereas the uncertain velocity changes are approximated using an sG particle projection.To ensure spectral accuracy in the random space, our particle method is designed to leverage the general structure of collision kernels discussed in [2].To this aim, we introduce a regularized kernel that offers several advantages, including the avoidance of iterative techniques commonly employed in other particle simulation methods [5].Being based on particle reconstruction, the method preserves the nonengativity of the solution and the main physical properties.To validate our method, we conduct several tests, these include classical benchmarks such as BKW, Trubnikov's solution for Coulombian particles, and the bump-on-tail problem.
The rest of the manuscript is organized as follows.In Section 2 we recall the basic ideas in the design of collision algorithm for the Landau-Fokker-Planck equation inspired by the corresponding grazing limit of the Boltzmann equation.Section 3 is devoted to present the details of the DSMC approach with a regularized kernel and the corresponding sG projection.A suite of numerical examples is then reported in Section 4 which confirms the efficiency and the accuracy of the present approach.The last section collects some concluding remarks and future developments.

Foundations of collision algorithms for the Landau equation
In this section, following [2,5], we recall some results concerning the derivation of direct simulation collision algorithm for the Landau-Fokker-Planck equation based on the grazing collision limit of the Boltzmann equation.To this aim, we will make use of a first order approximation of the space homogeneous Boltzmann equation and devote particular attention to the Maxwellian and the Coulombian cases.
We consider the space homogeneous Boltzmann equation [11] ∂f describing the time evolution of the one-particle distribution function f (v, t) with velocity v ∈ R 3 , at the time t > 0. The collision term Q(f, f ) on the right-hand side of (1) is a bilinear operator describing the binary interactions between particles and it is given by where q = v − v * is the relative velocity and n ∈ S 2 is the unit vector normal to the unit sphere S 2 in R 3 .The binary interactions rules characterizing collisions between particles are and the collisional kernel B(|q|, cos θ) reads where cos θ = q • n/|q| and σ(|q|, θ) is the so-called differential collision cross section at the scattering angle θ, corresponding to the number of particles scattered per unit of incident flux, per unit of solid angle, in the unit time.We introduce also the total scattering cross section and the momentum-transfer (or -transport) scattering cross section defined as that describes the average momentum transferred in the collisions.The differential cross section σ(|q|, θ) takes different forms depending on the interactions considered.Amongst the most relevant cases we mention the variable hard sphere (VHS) model and the Coulomb interactions.The VHS model is such that where C γ > 0 is a constant, with γ = 0 for Maxwell molecules and γ = 1 for hard sphere model, which gives On the other hand, particles subject to Coulomb forces are characterized, according to the Rutherford formula, by the scattering where e is the charge of the particles, ε 0 the vacuum permittivity and m r is the reduce mass, corresponding to m/2 for particles of the same species.In this case, it is necessary to introduce a cut-off, since for θ → 0 the cross section is singular.Following the usual approximations justified by the shielding effect, we have is the Debye length and It is well-known that in the grazing collision limit, from equation (1) with collisional operator (2), we can recover the Landau-Fokker-Planck equation [26] ∂f where Φ is a 3 × 3 nonnegative symmetric matrix and its usual form is given by Φ(q) = |q| γ+2 S(q), S(q) = I − q ⊗ q |q| 2 , with I identity matrix, and where γ = 0 is the Maxwellian case and γ = −3 is the Coulombian case.In the following, we will discuss in more detail this aspect and derive first order approximation formulas which are suitable for the design of DSMC type algorithms.

First order approximation of the Boltzmann equation
From a direct inspection of the collisional kernel in the presence of Coulomb forces, it is evident that when collisions become grazing, the interaction rate becomes infinite.To overcome this problem, which prevents the direct application of the DSMC algorithm, following [2] we introduce a suitable first order approximation of the Boltzmann collision operator that permits to avoid such restrictions.In the sequel, we omit the explicit dependence on the time variable for compactness of notation.First, we rewrite equation (1) as where U = (v + v * )/2 is the velocity of the centre-of-mass, and The operator J is defined such that its action on the angular variable ω = q/|q| is Now we expand the operator J at first order in ε as where Î is the identity operator and ε is a small parameter.Note that in the limit ε → 0 we recover the definition of J.The collision term is then approximated by where The action of the operator exp(εJ) on any test function ψ(ω) can be represented as the Green function, P l (ω • n) the Legendre polynomial of order l, and The approximated gain operator then reads Let us now assume that the differential scattering cross section σ(|q|, θ) is concentrated at a small angle θ ≈ 0, that is the grazing collision regime, thus µ in B ε (µ, |q|) is concentrated near the value µ = 1.We have where P l (1) = l(l + 1)/2.As a consequence, the Green function becomes Finally, one can show that equation (6) with the above collisional kernel approximates the Landau-Fokker-Planck equation for small values of ε.If we define and the approximated Green function with the cut-off of the scattering angle as we have the first order approximation where the first term on the right-hand side of the previous equation plays the role of the gain operator with the approximated Green function.
We consider now the cases of Maxwellian and Coulomb potentials in a single species plasma model.If γ = 0, then we have that that is independent from the relative velocity.If γ = −3, we recover In the following, we will refer to τ (γ) 0 simply as τ 0 , in the general case.

Approximated scattering kernels
In order to solve numerically equation (8), we need to sample the probability density D(µ, τ 0 ) defined in (7).To this end, a random sample of the quantity µ(ε) in [−1, 1], at each time step, is required.It is easily observed that this can be very heavy from a numerical perspective.However, as pointed out in [2], we may consider simpler collisional kernels D * (µ, τ 0 ) with the following properties 1. D * (µ, τ 0 ) ≥ 0, and 2π A general form of such a kernel is furnished by the following result.

Proof.
For the proof, we refer to [2].
In the following, we consider three different kernels.The first one has been proposed in [5,29] and is commonly used in plasma physics simulation, and it reads As a consequence, to sample D (1) * (µ, τ 0 ) in the spherical coordinate system (θ, φ), we need to solve at every time step the nonlinear equation ( 11) and then compute the scattering angle cos θ as cos θ and φ as where r 1 , r 2 are two uniform random numbers in (0, 1).
The second kernel has been proposed in [2] as a simplification of where Sampling D (2) * (µ, τ 0 ) corresponds, in the spherical coordinate, to compute φ as in ( 13) and to fix cos θ = ν(τ 0 ).(15) Note that in this case we need only one random number, i.e. r 2 , and no iterations are required.
A third approximation consistent with Lemma 1 will be considered in the following and it consists in a regularization of the previous kernel D (2) * (µ, τ 0 ) defined as We highlight that ν(τ 0 ) is obtained by fixing whereas φ is computed as in (13).It is easy to verify that the kernel D (3) * (µ, τ 0 ) satisfies conditions 1-2-3.

Particle methods and their sG projection
We introduce a random vector z = (z 1 , . . ., z d ) ∈ R dz collecting all the uncertainties of the system and for which it is known its distribution p(z) for any I z ∈ R dz .The vector z characterizes the uncertainties of the system due, e.g., to missing information on the initial state, measurement of the model parameters and, boundary conditions [27,46].
We are interested in the evolution of the density f (v, t, z), v ∈ R 3 , z ∈ R dz , t ≥ 0, solution to the Landau-Fokker-Planck equation ( 5) with uncertainties In particular, in (18) the interaction matrix Φ(•, z) may encode different types of interactions.
In this section, we first recall the basic ingredients of DSMC particle method in the absence of uncertainties, as proposed in [5].
Once established the DSMC solver of interest we seek to derive the related sG particle approach.These methods have been previously studied for mean-field model for emerging phenomena [9,10,28], and subsequently extended to the space-homogeneous Boltzmann equation in [36] and to a non-homogeneous plasma dynamics with BGK collision operator in [27].

DSMC methods
Let us now give the details of the DSMC algorithms based on the previous first order approximation of the Boltzmann equation in the grazing collision limit and in absence of uncertainties (see [2,5,29,32,45] and the references therein).
Hence, indicating with q n = v n i − v n j the relative velocity at the time step n of any two particles v n i , v n j , the collision law reads with h n given by where q n ⊥ = (q n y ) 2 + (q n z ) 2 1/2 .In the above formulas, the angles θ and φ are given by ( 12)-( 13), ( 15)-( 13), or ( 17)-( 13) according to the approximated kernel we consider.
We have denoted by Sround(x) the stochastic rounding of a positive real number x

Sround =
x + 1 with probability x − x x with probability 1 − x + x , where x denotes the integer part of x.In Algorithm 1 we recall the standard Nanbu-Babovsky algorithm for Maxwell particles with collisions defined by (19).

Bird's scheme
We can also construct a Monte Carlo scheme based on the classical Bird's no time counter method.In this case, the average number of collisions at each time step is given by so that the average time between two collisions is Note that, in our setting, Bird's method has the same structure of the Nanbu-Babovsky method where at each ∆t c only one pair collides, since As a consequence, recollision between particles are admissed in a time step ∆t = N ∆t c in contrast to Nanbu-Babovsky.In Algorithm 2 we recall the standard Bird's scheme for Maxwell particles with collisions defined by (19) in the absence of uncertainties.

Stochastic Galerkin reformulation of the DSMC method
We consider a sample of N particles v i (t n , z), i = 1, . . ., N at time t n = n∆t, and we expand them by their generalized polynomial chaos (gPC) expansion being {Ψ k (z)} M k=0 a set of polynomials of degree less or equal to M ∈ N, orthonormal with respect to the measure p(z)dz where Ω ∈ R d and δ kl is the Kronecker delta function.The polynomials {Ψ h (z)} M k=0 are chosen following the so-called Wiener-Askey scheme, depending on the distribution of the parameters.In (20), for fixed k = 0, . . ., M , the term vn i,h is the projection in the space generated by the polynomial of degree k ≥ 0 We can also write where we indicated with qn ij,k (t) the following quantity We recall that in general τ = τ (z), since it may depend on the relative velocity, and then we also have τ 0 = τ 0 (z) and θ = θ ij (z).We substitute the approximation of the velocities into (19) and we project against Ψ l (z)p(z)dz for every l = 0, . . ., M to obtain where the collision matrices are defined as Therefore, we may rephrase the Nanbu-Babovsky scheme in the sG framework as presented in Algorithm 3. Whereas the sG version of the Bird's scheme is given in Algorithm 4.

Algorithm 3 sG Nanbu-Babovky for the Landau equation
• Compute the initial gPC expansions {v M,0 i } N i=1 from the initial distribution f 0 (v); • for n = 1 to nTOT, given the projections {v n i,k , i = 1, . . ., N, k = 0, . . ., M }: set Nc = Sround(ρN ∆t/2τ ); select the interaction pairs (i, j) uniformly among all the possible ones, and for every pair (v n i,k , vn j,k ): * compute the matrices Vkl and Ŵk according to (23); * perform the collision according to (22); * set vn+1 i,k = v i,k and vn+1 j,k = v j,k ; set vn+1 i,k = vn i,k for all the particles that have not been collided; • end for.

Numerical examples and applications
In this section, we present several numerical test and examples to validate our algorithms both without and with uncertain parameters.First, we investigate the Maxwellian case without uncertainties.Then, we show several tests for the DSMC-sG method, for both the Nanbu-Babovsky and the Bird's schemes.We check for the spectral convergence and the accordance with the exact BKW solution of the Maxwellian case with uncertainties.Then, we consider the Coulombian case, focusing the attention on the regularity of the kernels D (1) * , D (2) * , and D (3) * .We check again for the spectral convergence, and then we test the schemes on the standard case studies for the homogeneous Landau equation.In particular, we concentrate on the capability to reach the equilibrium starting from different uncertain initial conditions, i.e., anisotropic initial temperature, sum of two Gaussians, and bump-ontail distributions.In all the subsequent tests, we consider the physical constants fixed as follows: e = ε 0 = m = ρ = 1, and log Λ = 0.5.

Test 1: Exact solution in the Maxwellian case
We consider the model with Maxwell molecules, i.e., γ = 0 in (3).In 3D, an exact solution is given by (see [7], Appendix A) with T temperature and We recall that the moments of order zero, one, and two are conserved, while the exact time evolution of the fourth order moment reads In Figure 1, we show the relative L 2 errors of the distribution f (v, t) and the fourth order moment M4(t) of the DSMC approximation with respect to the exact BKW solution (24)   27) and ( 28).The particles are N = 5×10 6 , the time step ∆t = 0.1, and the temperature T = 0.07.We choose the Nanbu-Babovsky scheme with the kernel D (3) * .Initial conditions given by ( 26), with z-temperature T 0 z = 0.04 for the Maxwellian case, and T 0 z = 0.04, 0.001 for the Coulombian case.
In Figure 2 we display the time evolution of the first and second order moments (left panel), and the fourth order moment (right panel) computed with the DSMC scheme, together with the exact solutions.We choose the Nanbu-Babovsky algorithm, kernel D (3) * , N = 5 × 10 6 particles, ∆t = ε/ρ = 0.1, and initial conditions given by (24) with t = 0 and T = 1.

Test 2: Trubnikov test
We initialize the distribution as an ellipsoid, i.e., as with anisotropic initial temperature The temperature difference ∆T (t) = T x (t) − T z (t) goes to zero exponentially with a specific rate τ T ∆T (t) = ∆T (0)e −t/τ T .
In the Maxwellian case (see Appendix A for further details), we have in the Coulombian case Trubnikov [41] obtained an approximated solution in the limit of small temperature difference |T 0 We choose N = 5 × 10 6 particles, ∆t = 0.1, and the Nanbu-Babovsky scheme with the kernel D (3) * .We fix the total temperature T = 0.07 and we vary the initial z-temperature T 0 z to change ∆T (0).In Figure 3 we compare the benchmark solutions with the particle approximations, for both the Maxwellian and the Coulombian cases.In the first scenario (left panel), the decreasing rate does not depend on the magnitude of the temperature difference, therefore we choose T 0 z = 0.04 so that ∆T (0) = 0.045 but any other choice gives the same result.In the Coulombian case (right panel), we choose T 0 z = 0.04, 0.001 to have ∆T (0) = 0.045, 0.1035.As we can see, for higher ∆T (0) the decreasing rate is not in accordance with the Trubnikov solution.

Test 3: Maxwellian case with uncertainties
We consider here the model with Maxwell molecules and uncertainties in the initial temperature T (z), i.e., with and fourth order moment given by We initialize the particles in the DSMC-sG scheme by sampling the initial conditions given by   In all the tests, we adopt the Nanbu-Babovsky scheme given by Algorithm 3 and the kernel D In Figure 4 we show the stochastic Galerkin error for increasing values of M (left) and its time evolution at fixed M (right).We compute a reference solution with M = 30 and we store the collisional tree.Then, for different values of M , we compute the L 2 -Error of the temperature T (z) with respect to the reference solution.We notice that the machine accuracy is reached with a finite number of modes, and the error is constant in time.
In Figures 5-6-7 we may see the comparison between the BKW exact solution given by ( 29) and the DSMC-sG approximation.We choose N = 5 × 10 7 particles and a stochastic Galerkin expansion up to order M = 5.In particular, in Figure 5 we display at fixed times t = 0, 1, 5 the marginals of E z [f (v, t, z)] and Var z [f (v, t, z)].We may notice the accordance between the DSMC-sG approximation (red circles) and the exact BKW solution (solid black lines).In Figures 6-7 the show at the fixed times t = 0, 1, 5 the slices of the three dimensional plots of E z [f (v, t, z)] and Var z [f (v, t, z)] respectively.In both the figures, the top row is the DSMC-sG approximation, while the bottom row is the exact BKW solution.We observe again a good agreement.

Stochastic Galerkin error
We are interested first in evaluating the convergence of the scheme in the space of the random parameters.To this aim, we consider the initial conditions f 0 (v, z) = f (v, 0, z) where the total temperature, conserved in time, is defined as with uncertain initial conditions in the temperature along the x axis and z ∼ U([0, 1]).We choose ∆t = 0.1, in a way that ε = ρ∆t = 0.1, and N = 10 6 particles.We consider a reference solution obtained with M = 30 up to time t = 1 and we store both the initial data and the collisional tree.Then, we compute the L 2 error in the evaluation of the total temperature T (z) for increasing M .The results obtained with the Nanbu-Babovky (NB) and the Bird (B) scheme are similar, therefore we show only the NB results.
In Figure 8 we show the decay of the L 2 error in the evaluation of the temperature T (z), for the kernels D (3) * .We observe that the error convergence deteriorates for D (1) * , since we need to introduce a cut-off for the numerical resolution of the nonlinear equation (11), and for D (2) * , because the function ν(τ 0 ) is not differentiable for τ 0 = 1.With the choice D (3) * , we reach the machine accuracy with a small order M .

Trubnikov test
We initialize the distribution as an ellipsoid, i.e., as , with anisotropic uncertain initial temperature (right panels).The black line is the benchmark solution (32).The number of particles is N = 10 6 , M = 5 and ∆t = ε/ρ.
In the case of small uncertain temperature difference, the Trubnikov [41] formula reads ∆T (z, t) = ∆T 0 (z) exp (−t/τ T ) , with τ T = 5 8 which can be used as a benchmark solution to test the numerical schemes.We choose T 0 x (z) = T 0 y (z) = 0.08 + 0.04z T 0 z = 0.04, with z ∼ U([0, 1]), N = 10 6 particles and a sG expansion up to order M = 5.First, we analyse the behaviour of the schemes for different values of the scale parameter ε = ρ∆t.We test the Nanbu-Babovsky and the Bird scheme, for the kernels D (3) * .The results are summarised in Figure 9.We observe that for smaller values of ε, the numerical solutions get closer to the reference one, indicating that the schemes are consistent with the grazing limit.However, as pointed out also in [5], there is a discrepancy between the two solutions, since the Trubnikov relation ( 32) is still an (analytical) approximation.
Then, at fixed ε = 0.5, we compare the solutions obtained with different kernels, for both the Nanbu-Babovsky and the Nanbu scheme.In Figure 10 we may notice that with the choices D In the end, we compare the results obtained with the Nanbu-Babovsky and the Nanbu schemes for all the kernels.We fix ε = 0.5, N = 10 6 particles and M = 5.In Figure 11 we show the time evolution of the temperature difference for D

Sum of two Gaussians
We initialize the distribution as the sum of two Gaussians centred in v = ±1 and with the same uncertain temperature T (z) = 0.1 + 0.2z, z ∼ U([0, 1]), * , and D (3) * , for fixed ε = 0.5, for both the Nanbu-Babovsky (left panel) and Bird (right panel) schemes.The black line is the benchmark solution (32).The number of particles is N = 10 6 , M = 5 and ∆t = ε/ρ = 1.i.e., We use the Nanbu-Babovsky algorithm and the kernel D (3) * .We choose N = 5 × 10 7 particles, ∆t = ε/ρ = 1 and a stochastic Galerkin expansion with M = 5.As we expected, the system evolves toward the equilibrium distribution, that is the centred Gaussian with the temperature T (z).In Figure 12 we show the slices at fixed times t = 0, 100, 1000 of E z [f (v, t, z)] (upper row), and Var z [f (v, t, z)] (lower row).In Figure 13 we display the marginals of E z [f (v, t, z)] at the same times.From the right panel, we may notice the accordance between the numerical solution (black circled line) and the expected equilibrium distribution (red starred line).

Bump on tail
We initialize the distribution as a centred Gaussian with a bump on tail, namely a small portion of mass concentrated in v = 3 with uncertain temperatures T 1 (z) = 0.2 + 0.2z, z ∼ U([0, 1]), and T 2 (z) = T 1 (z)/40.Given this initial condition, we observe that the conserved quantities, besides the mass, are M1 =  In Figure 14 we display the marginals E z [f (v, t, z)] at fixed times t = 0, 20, 50, 500.We observe that the bump is absorbed as the time increase and, at the time t = 500, the system is close to the equilibrium, i.e. the Gaussian with mean M1 = 3  40 and temperature T (z).In particular, in the bottom row, right panel, the accordance is good between the DSMC-sG approximation and the expected equilibrium distribution.

Test 5: DSMC-sG versus DSMC-MC
We are interested here in comparing the DSMC-sG algorithm with the DSMC-MC, i.e., the standard DSMC scheme with a Monte Carlo sampling of the random parameter.We take the same computational setting of Section 4.3, that is, Nanbu-Babovsky scheme, kernel D (3) * , and initial conditions given by (31).We fix the number of particles N = 10 6 , the time step ∆t = ε/ρ = 0.1, and we compute the error in the evaluation of the fourth order moment at fixed time t where M4(t, z) is a reference solution computed with a collocation algorithm with N = 10 8 particles and M = 20 collocation nodes.We want to compare the error with the computational cost, that is N • M 2 for the DSMC-sG, where M is the order of expansion, and N • M for the DSMC-MC, where M is the size of the sample of the random parameter.In Figure 15 we show the results for the cost divided by the fixed number of particles N .As we can see, the DSMC-sG performs better than the DSMC-MC for small cost.Moreover, we observe that the error of the DSMC-sG is saturated by the Monte Carlo part of the algorithm since it is constant for increasing order of expansion. -

Conclusions
In this work, we have investigated the design of efficient particle methods for the Landau-Fokker-Planck equation in the presence of uncertainties.The equation has significant implications for the development of fusion reactors and has been extensively utilized by researchers and engineers to gain insights into the behaviour of charged particles in plasma.To accurately predict the behaviour of these particles, it is crucial to account for uncertainties in the constitutive parameters that characterize the system behaviour.
Our approach combines collision algorithms inspired by the grazing limit of the Boltzmann equation with a stochastic Galerkin particle projection.This method leverages the general structure of kernels in collision algorithms and introduces a regularized kernel to ensure spectral accuracy the random space.In addition, the method benefits from a more efficient collision strategy which avoids iterative algorithms and, by employing particle reconstruction, it preserves the nonnegativity of the solution and other essential physical properties.
To validate the effectiveness and accuracy of our approach, we have conducted various numerical tests, including classical benchmarks such as BKW, Trubnikov's solution for Coulombian particles, and the bump-on-tail problem.The results demonstrate the efficiency and accuracy of our method in capturing the complex dynamics described by the Landau equation with uncertain data.
In conclusion, we have developed a robust and accurate approach that can effectively capture the behaviour of charged particles in collisional plasmas under uncertain data.Future developments may involve extending the methodology to more complex scenarios and investigating additional applications in plasma physics in space non-homogeneous situations.

A Trubnikov formula for Maxwell molecules
We follow [41] (Section 20, pages 200-203) to derive the Trubnikov formula for the relaxation of anisotropic initial temperatures in a system characterized by Maxwell molecules.We consider equation ( 5) with γ = 0 written in flux form with initial conditions given by where T ⊥ = T x = T y > T z and v 2 ⊥ = v 2 x + v 2 y .Since the temperature T = (2T ⊥ + T z )/3 is constant in time, we have where J z is the z-component of the flux Using the expression of the initial distribution, we can compute explicitly the sum inside the integral, which gives In the Coulombian case, we have a multiplicative 1/q 3 term inside the double integral, which forces us to make the assumption that |T ⊥ − T z | 1, (see [41], Section 20, page 202, for further details).In the Maxwell case, we can compute the double integral explicitly, using again the expression of the initial distribution.Through the change of variables (v, v * ) → (V, q), where V = (v + v * )/2 and q = v − v * , we have Plugging the previous expression into (35) we get We remark that this expression holds independently of the magnitude of |T ⊥ − T z |.Furthermore, the rate inside the exponential function does not depend on the temperature T as in the Coulombian case.

Figure 1 :
Figure 1: Test 1.Time evolution of the relative L 2 errors with respect to the BKW solution of the distribution f (v, t) (left) and of the fourth order moment M4(t) (right), for the kernels D (1) * , D (2) * , and D (3) * , and different values of ε = ρ∆t.In all the tests, we use 5 × 10 6 particles and the Nanbu-Babovsky scheme.Initial conditions given by (24) with t = 0 and T = 1.

Figure 2 : 1 .Figure 3 :
Figure 2: Test 1.Comparison between the time evolution of the first and second order moments (left), and of the fourth order moment (right), obtained with the DSMC method and from the exact BKW solution.In all the tests, we use 5 × 10 6 particles, ∆t = ε/ρ = 0.1, the kernel D (3) * , and the Nanbu-Babovsky scheme.Initial conditions given by (24) with t = 0 and T = 1.

6 Figure 4 :
Figure 4: Test 3. Left: L 2 -Error in the evaluation of the temperature T (z) at fixed time t = 1 forincreasing M , with respect to a reference solution, for different values of κ.Right: time evolution of the same error in the time span [0, 5] for the case κ = 0.95.We consider in both cases N = 10 6 particles, ∆t = ε/ρ = 0.1, and initial conditions given by(31).The kernel is D(3) * , the scheme is Nanbu-Babovsky.Reference solution computed with M = 30.

Figure 8 :
Figure 8: Test 4. Comparison of the L 2 error in the evaluation of the temperature T (z) for different kernels, for increasing M , with the NB scheme.We choose N = 10 6 , t = 1, ∆t = ε/ρ = 0.1 and the reference solution is computed with an order M = 30.

( 1 )
* , D(2) * , and D better with respect to the choice D (1) * , in the sense that they are closer to the Trubnikov solution.

Figure 10 :
Figure 10: Test 4. Comparison of the time evolution of the expectation of ∆T (z, t)/∆T (z, 0) for the choices of the kernel D (1) * , D

Figure 13 :
Figure 13: Test 4. Evolution at fixed times t = 0, 100, 1000 of the marginal E z [f (v, t, z)] of the Coulombian model with uncertainties, with sum of two Gaussians initial conditions (33).We choose N = 5 × 10 7 particles, M = 5, ∆t = ε/ρ = 1.The equilibrium marginal distribution (red starred) of the right panel is the Maxwellian distribution computed with the theoretical (conserved) mean and temperature.

Figure 14 :Figure 15 :
Figure 14: Test 4. Evolution at fixed times t = 0, 20, 50, 500 of the marginal E z [f (v, t, z)] of the Coulombian model with uncertainties, with bump on tail initial conditions.We choose N = 5 × 10 7 particles, M = 5, ∆t = ε/ρ = 1.The equilibrium marginal distribution (red starred) of the lower row, the right panel is Maxwellian distribution computed with the theoretical (conserved) mean and temperature.