Stochastic variational principles for the collisional Vlasov-Maxwell and Vlasov-Poisson equations.

In this work, we recast the collisional Vlasov-Maxwell and Vlasov-Poisson equations as systems of coupled stochastic and partial differential equations, and we derive stochastic variational principles which underlie such reformulations. We also propose a stochastic particle method for the collisional Vlasov-Maxwell equations and provide a variational characterization of it, which can be used as a basis for a further development of stochastic structure-preserving particle-in-cell integrators.

Introduction

The collisional Vlasov equation describes the time evolution of the particle density function of plasma consisting of charged particles of charge q and mass m which undergo collisions described by the collision operator , and are subject to the electric and magnetic fields. The vectors and denote positions and velocities, respectively. For simplicity, we restrict ourselves to one-piece plasmas. Usually, the particle density function is normalized, so that the total number of particles is . However, in this work we would like to treat f as a probability density function, and therefore we will use the normalization instead. A self-consistent model of plasma is obtained by coupling (1.1) with the Maxwell equations where denote the charge density and the electric current density, respectively, and the factor is due to our normalization. The system (1.1)–(1.3) is usually referred to as the Vlasov–Maxwell equations. It will also be convenient to express the electric and magnetic fields in terms of the scalar and vector potentials and as is typical in electrodynamics. The Vlasov–Poisson equations are an approximation of the Vlasov–Maxwell equations in the non-relativistic zero-magnetic field limit (see §6). The main goal of this work is to provide a variational characterization of the Vlasov–Maxwell and Vlasov–Poisson equations via a stochastic Lagrange–d’Alembert type of a principle.

Variational principles have proved extremely useful in the study of nonlinear evolution partial differential equations (PDEs). For instance, they often provide physical insights into the problem being considered; facilitate discovery of conserved quantities by relating them to symmetries via Noether’s theorem; allow one to determine approximate solutions to PDEs by minimizing the action functional over a class of test functions (e.g. [1]); and provide a way to construct a class of numerical methods called variational integrators [2,3]. A variational principle for the collisionless Vlasov–Maxwell equations was first proposed in [4]. It has been used to derive various particle discretizations of the Vlasov–Maxwell and Vlasov–Poisson equations [5-9], including structure-preserving variational particle-in-cell (PIC) methods [10-12]. It has also been applied to gyrokinetic theory (e.g. [13,14]). For other formulations and extensions, see also [15].

A structure-preserving description of collisional effects is far less developed. A metriplectic framework for the Vlasov–Maxwell-Landau equations has been presented in [16,17]. More recently, a stochastic variational principle has been proposed in [18] to describe collisional effects for the Vlasov equation with a fixed external electric field. To the best of our knowledge, to date no variational principle has been derived for the collisional Vlasov–Maxwell and Vlasov–Poisson equations. In this work, we extend the notion of the stochastic Lagrange–d’Alembert principle presented in [18] to plasmas evolving in self-consistent electromagnetic fields. The main idea of our approach is to interpret the Vlasov equation (1.1) as a Fokker–Planck equation and consider the associated stochastic differential equations.

The idea of using stochastic differential equations to model collisions has been pursued by a number of authors over the last few decades (e.g. [18-32], Y Fu, X Zhang, H Qin 2020, unpublished data).

There has been an ever-growing body of the literature dedicated to stochastic variational principles in recent years. Stochastic variational principles allow the introduction of noise into systems in such a way that the resulting probabilistic models retain all or some of the geometric properties of their deterministic counterparts. For this reason, stochastic variational principles have been considered in the context of Lagrangian and Hamiltonian mechanics [18,33-39], soliton dynamics [40,41], fluid dynamics [42-49] and kinetic plasma theory [18].

Main content. The main content of the remainder of this paper is, as follows.

In §2, we recast the collisional Vlasov–Maxwell equations as a system of coupled stochastic and partial differential equations.

In §3, we discuss the relationship between particle methods and stochastic modelling. We formulate a stochastic particle discretization for the collisional Vlasov–Maxwell equations and cast it in a form that allows the derivation of a variational principle.

In §4, we describe the variational structure underlying the stochastic particle discretization of the Vlasov–Maxwell system. The main result of this section is theorem 4.2, in which a stochastic Lagrange–d’Alembert principle for the particle discretization is proved.

In §5, we generalize the ideas from §4 to the original undiscretized equations. The main result of this section is theorem 5.1, in which a stochastic Lagrange–d’Alembert principle is proved for a class of the collisional Vlasov–Maxwell equations.

In §6, we prove a stochastic Lagrange–d’Alembert principle applicable to the Vlasov–Poisson equations. The main result of this section is theorem 6.1.

Section 7 contains the summary of our work.

The Vlasov–Maxwell–Fokker–Planck equations

Stochastic reformulation

Various collision models and various forms of the collision operator are considered in the plasma physics literature (e.g. [50,51]). A key step towards a stochastic variational principle is a probabilistic interpretation of the Vlasov equation (1.1). Therefore, in this work we will be interested only in those collision operators for which (1.1) takes the form of a linear or strongly nonlinear Fokker–Planck equation (e.g. [52-54]). Namely, we will assume that the collision operator can be expressed as for some symmetric positive semi-definite matrix and vector functions, where the dependence of and on f may in general be nonlinear, and may involve differential and integral forms of f. In that case (1.1) is an integro-differential equation, the so-called strongly nonlinear Fokker–Planck equation [52]. In case and are independent of f, that is, and , the Vlasov equation (1.1) reduces to the standard linear Fokker–Planck equation. We will further assume that and can be expressed in the form for a vector function , and a family of vector functions with . Note that given a symmetric positive semi-definite matrix , a decomposition (2.2) can always be found, but it may not be unique. For instance, one may take and assume that for . Then the first equation in (2.2) implies that the family of functions can be determined by calculating the square root of the matrix , and the second equation in (2.2) can be used to calculate the function G. If (1.1) has the form of a Fokker–Planck equation, then the particle density function f can be interpreted as the probability density function for a stochastic process . This stochastic process then satisfies the Stratonovich stochastic differential equation [52-55] and where denote the components of the standard M-dimensional Wiener process and ○ denotes Stratonovich integration. Note that the terms G and can be interpreted as external forces, and that in their absence the equations (2.3) reduce to the equations of motion of a charged particle in an electromagnetic field. We will therefore refer to G and as forcing terms. The electric and magnetic fields are coupled via the Maxwell equations (1.2). It should also be noted that unless (1.1) is linear, the right-hand side of (2.3) depends on f. In order to obtain a self-consistent system, one can express f in terms of the stochastic processes X and V as , where denotes the expected value, and is Dirac’s delta. This can be further plugged into (1.3). Together, we get Equations (1.2), (2.3) and (2.4) form a self-consistent system of stochastic and partial differential equations whose solutions are the stochastic processes , , and the functions , .

Remark.

Upon substituting (2.4a), the forcing terms G and become functionals of the processes X and V, that is, and . However, for convenience and simplicity, throughout this work we will stick to the notation and , understanding that the probability density is given by (2.4a) (or by (3.2a) for particle discretizations; see §3).

Examples

Below we list a few examples of collision operators that fit the description presented in §2a.

Lenard–Bernstein operator

The Lenard–Bernstein collision operator where , and are parameters, models small-angle collisions, and was originally used to study longitudinal plasma oscillations [50,51,56]. It can be easily verified that an example decomposition (2.2) for is given by the functions Note that these functions do not explicitly depend on f, therefore in this case (1.1) is a linear Fokker–Planck equation.

Lorentz operator

The Lorentz collision operator models electron–ion interactions via pitch-angle scattering and is given by the formula where is the collisional frequency as a function of the absolute value of velocity, is the identity matrix and denotes tensor product. The primary effect of this type of scattering is a change of the direction of the electron’s velocity with negligible energy loss. More information about the Lorentz collision operator, including the exact form of the collision frequency, can be found in, e.g. [50,51,57,58]. It can be verified by a straightforward calculation that an example decomposition (2.2) for is given by the functions Note that these functions do not explicitly depend on f, therefore also in this case (1.1) is a linear Fokker–Planck equation.

Coulomb/Landau operator

The more general Coulomb collision operator has the form (2.1) with where appears due to our normalization of f, is Kronecker’s delta, and , with denoting the so-called Coulomb logarithm. The Coulomb operator describes collisions in which the fundamental two-body force obeys an inverse square law, and makes the assumption that small-angle collisions are more important than collisions resulting in large momentum changes [50,51,59]. A decomposition (2.2) can be found, for example, via the procedure outlined in §2a. However, the expressions for G and are complicated, therefore we are not stating them here explicitly. Note that and explicitly depend on f. Therefore, for the Coulomb operator the Vlasov equation (1.1) is a strongly nonlinear Fokker–Planck equation. Note also that and can be explicitly written as functionals of the stochastic processes X and V as The collision operator (2.1) with and as in (2.9) can also be expressed in an equivalent, although more symmetric form, known as the Landau form of the Coulomb operator, or simply the Landau collision operator (e.g. [50]).

Stochastic particle discretization of the Vlasov–Maxwell equations

Particle modelling is one of the most popular numerical techniques for solving the Vlasov equation (e.g. [60,61]). In this section, we discuss the connections between particle methods and stochastic modelling.

The standard particle method for the collisionless Vlasov equation (1.1) (with ) consists of substituting the Ansatz for the particle density function, and deriving the corresponding ordinary differential equations satisfied by the ‘particle’ positions and velocities , which turn out to be the characteristic equations. Note that we did a qualitatively similar thing in §2a, where we turned the original collisional Vlasov equation into the system of stochastic differential equations (2.3), which in the absence of the forcing terms G and have the same form as the characteristic equations, and in fact the ‘particles’ and can be interpreted as realizations of the stochastic processes and for different elementary events .

When the right-hand side of (2.3) does not depend on f, then (2.3) can in principle be solved numerically with the help of any standard stochastic numerical method (e.g. [55]), and each realization of the stochastic processes can be simulated independently of others. When the right-hand side of (2.3) depends on f, then all realizations of the stochastic processes have to be solved for simultaneously, so that at each time step the probability density function f can be numerically approximated (e.g. [52]). Such an approach, however, does not quite lend itself to a geometric formulation. Therefore, in order to be able to introduce a variational principle in §4, let us consider stochastic processes , with each pair satisfying the stochastic differential system and for , where are N independent M-dimensional Wiener processes. Note that the systems (3.1) are decoupled from each other for different values of , and each system is driven by an independent Wiener process . Therefore, the pairs for are independent identically distributed stochastic processes, each with the probability density function f that satisfies the original Fokker–Planck equation (1.1). In that sense (3.1) is equivalent to (2.3). The advantage is that instead of considering N realizations of the six-dimensional stochastic process in (2.3), one can consider one realization of the -dimensional process in (3.1). Such a reformulation will allow us to identify an underlying stochastic variational principle in §4. The last step leading to the stochastic particle discretization is approximating the probability density function f in (3.1). This can be done with the help of the law of large numbers, namely, one can approximate (2.4) for large N as It is easy to see that (3.2a) coincides with the standard Ansatz used in particle modelling (with the weights ). Therefore, the system of stochastic differential equations (3.1) with the approximation (3.2), and with the electromagnetic field coupled via the Maxwell equations (1.2), can be considered as a stochastic particle discretization of the collisional Vlasov–Maxwell equations.

Remark.

Upon substituting (3.2a), the forcing terms G and become functionals of the processes and . Similar to the discussion in §2a, for convenience and simplicity, throughout this work we will stick to the notation and , understanding that the probability density is given by (3.2a) for particle discretizations.

Variational principle for the particle discretization

In this section, we propose an action functional which can be understood as a stochastic version of the Low action functional [4], and we prove a variational principle underlying the particle discretization introduced in §3, akin to the stochastic Lagrange–d’Alembert principle first introduced in [18].

Function spaces

Before we introduce the action functional, we need to identify suitable function spaces on which it will be defined. For simplicity, let our spatial domain be the whole three-dimensional space , and let us consider the time interval for some . Let be the probability space with the filtration , and let for denote N independent M-dimensional Wiener processes on that probability space (such that is -measurable for all ). The stochastic processes and satisfy (3.1), so they are in particular -adapted semimartingales, and have almost surely continuous paths [62]. We also notice that there is no diffusion term in (3.1), therefore we even have that the processes are almost surely of class . We introduce the notation Note that this set is a vector space [62]. The potentials and A satisfy the Maxwell equations (1.2) and (1.4), therefore we require them to be of class . However, since our spatial domain is unbounded, we further need to assume that the vector fields E and B are square integrable. We introduce the notation where is simply the space of compactly supported elements of .

Action functional

Let us consider the action functional defined by the formula where denotes the time derivative of , and the electric and magnetic fields E and B are expressed in terms of the partial derivatives of the potentials and A as in (1.4). Following the standard convention in stochastic analysis, we will omit writing elementary events as arguments of stochastic processes unless otherwise needed, i.e. . The action functional (4.4) resembles the Low action functional introduced in [4]. In fact, it can be viewed as a particle discretization of the Low action functional, written in terms of stochastic processes [5,6,8,10-12]. The term is the so-called Hamilton–Pontryagin kinematic constraint (e.g. [63,64]) that enforces that using the Lagrange multiplier , which turns out to be the conjugate momentum. In principle, this constraint is not necessary in our context—we could omit it and replace with in (4.4). We will, however, keep it in order to make a clear connection with the theory developed in [35]. It also makes the notation in the proof of the stochastic Lagrange–d’Alembert principle in §4c more convenient and elegant. Note that the action functional S is itself a random variable, as is one of its arguments. The variations of S with respect to its arguments are given by (see appendix A for the details of the derivations) where ρ and J are defined in (3.2b) and (3.2c), respectively. The total variation of S with respect to the variations of all arguments equals

The stochastic Lagrange–d’Alembert principle

While the standard rules of the calculus of variations apply to the variations (4.5d,e), the variations (4.5a–c) involve stochastic processes and stochastic integrals. Therefore, before we can formulate a stochastic variational principle, we need the following lemma, whose proof is given in appendix B.

Lemma 4.1.

Let and , and let be of class for . Then if and only if where ‘a.s.’ means almost surely.

Equation (4.8) means that , and satisfy a stochastic differential equation, which can be written in the differential form as

We are now in a position to formulate and prove a stochastic variational principle that generalizes the deterministic Lagrange–d’Alembert principle for forced Lagrangian and Hamiltonian systems, akin to the stochastic variational principle introduced in [18].

Theorem 4.2 (Stochastic Lagrange–d’Alembert principle for particles).

Let and for be stochastic processes, and let , be functions. Assume that and for are functions of their arguments, where f is given by (3.2a). Then , , , and satisfy the system of stochastic differential equations for and , together with the Maxwell equations (1.2), (1.4) and (3.2) on the time interval , if and only if they satisfy the following variational principle for arbitrary variations , , , and , with almost surely, and for all , where the action functional S is given by (4.4).

Proof.

Let us first consider the variations with respect to A in (4.11). Given the boundary conditions for , from the standard calculus of variations we have that (see equation (4.5d)) for all if and only if (1.2d) is satisfied. Similarly, (see equation (4.5e)) holds for all if and only if (1.2a) holds. Further, for variations with respect to we have that (see equation (4.5b)) for all if and only if (4.10b) is satisfied almost surely, which follows from the standard theorem of the calculus of variations, since the integral in (4.5b) is a standard Lebesgue integral, and the integrands are almost surely continuous. Similarly, (see equation (4.5c)) for all if and only if (4.10a) is satisfied almost surely. Finally, for variations with respect to , equations (4.5a) and (4.11) give which, by lemma 4.1, holds for all if and only if (4.10c) is satisfied.

Equation (4.10) is expressed in terms of the Lagrange multipliers , which, as can be seen in (4.10b), turn out to be the conjugate momenta. The conjugate momenta can be eliminated, and equation (4.10) can be recast as equation (3.1b), which is shown in the following theorem.

Theorem 4.3.

Equations (3.1) and (4.10) are equivalent.

By calculating the stochastic differential on both sides of (4.10b) and substituting (4.10a), we obtain for each and . Comparing this with (4.10c), and using (1.4), one eliminates the conjugate momenta and obtains equation (3.1d).

Theorems 4.2 and 4.3 provide a variational formulation of the stochastic particle method from §3. One can further perform a variational discretization of the electromagnetic fields A and , for instance along the lines of [10,65] or [66], thus obtaining a stochastic PIC discretization of the collisional Vlasov–Maxwell equations. The resulting structure-preserving numerical methods will be investigated in a follow-up work.

Variational principle for the Vlasov–Maxwell equations

The form of the action functional (4.4) and of the Lagrange–d’Alembert principle (4.11) suggests that it should be possible to formulate a similar variational principle for the stochastic reformulation of the Vlasov–Maxwell system discussed in §2a. In this section, we provide such a variational principle for a class of collision operators.

Let us consider the action functional defined by the formula where denotes the time derivative of X, the electric and magnetic fields E and B are expressed in terms of the partial derivatives of the potentials and A as in (1.4), and denotes the expected value of the random variable . Note that unlike S in (4.4), the action functional is not a random variable, as the dependence on is integrated out with respect to the probability measure by calculating the expected value. In fact, S could be regarded as a Monte Carlo approximation of when the processes are independent and identically distributed as X, and similarly for V and P. An important issue to consider is the domain of this action functional. In a similar manner to (4.3), one may want to take as the domain the set on which the formula (5.1) is well defined. This domain, however, turns out to be too big, in the sense that, as will be discussed below, due to the presence of the expected value the variations of do not uniquely determine the set of stochastic evolution equations (2.3). It is therefore necessary to restrict (5.2) to a smaller subspace or submanifold which is compatible with the considered collision operator. Below we will demonstrate how this can be done for a class of collision operators (2.1) for which , that is, we have This class encompasses, for instance, the Lenard–Bernstein operator (2.6), or the more general nonlinear energy and momentum preserving Dougherty collision operator and its modifications [67-74]. For a given collision operator of the form (5.3), we define a compatible subset of , namely, Note that for any we have that , that is, . Therefore, the pair is an affine subspace of . The action functional can now be defined as Similar to the calculations in §4b, the variations of with respect to V and P are given by, respectively, and except that here , so that . For the variation of with respect to X we have Since , we have that . Furthermore, the variations are almost surely of class , and therefore have sample paths of almost surely finite variation. Consequently, the quadratic covariation almost surely [62]. Since the expected value of the Itô integral with respect to the Wiener process is zero, we altogether have that By plugging this in (5.8), we finally obtain The variations with respect to A and are the same as in (4.5d,e), respectively, only with the charge and electric current densities given by (2.4) rather than (3.2). The total variation of with respect to the variations of all arguments is given by

In the following theorems, we establish a variational principle for the system of equations (1.2), (2.3) and (2.4) for a class of collision operators with for all .

Theorem 5.1 (Stochastic Lagrange–d’Alembert principle for the VM equations).

Let , , be stochastic processes, and let , be functions. Assume that is a function of its arguments, where f is given by (2.4a). Then , , , and satisfy the system of stochastic differential equations for , together with the Maxwell equations (1.2), (1.4) and (2.4) on the time interval , if and only if they satisfy the following variational principle for arbitrary variations , , , and , with almost surely, and for all , where the action functional is given by (5.1) and (5.5).

Similar to the proof of theorem 4.2, the equations and are equivalent to (1.2a) and (1.2d), respectively. Note that is a subspace of , and is an inner product on that space. Therefore, by substituting equations (5.6), (5.7) and (5.10) in equation (5.13), and using the fact that the variations are arbitrary, we establish equivalence with equations (5.12a,), as well as with the equation for , which in turn is equivalent to equation (5.12c), given the assumption .

Theorem 5.2.

Equation (2.3) with for and equation (5.12) are equivalent.

Similar to the proof of theorem 4.3, by calculating the stochastic differential on both sides of equation (5.12b) and comparing with equation (5.12c), one eliminates P and obtains equation (2.3b).

Note that the forcing terms do not explicitly appear in the variational equation (5.13). By comparing theorem 4.2 and theorem 5.1, one could intuitively expect that the relevant variational principle should read However, due to the presence of the expected value in this equation, part or all of the information about the Stratonovich integral term is lost, as we saw in (5.9) for instance. Therefore, if the domain (5.2) is chosen for , then the variational equations (5.13) or (5.15) do not determine a unique set of stochastic differential equations that need to be satisfied by the considered stochastic processes. Consequently, it is necessary to encode the missing information about the forcing terms in the definition of the action functional by restricting its domain to a subset compatible with the considered collision operator. For the class of collision operators (5.3) a suitable choice of the domain is proposed in (5.5). For other collision operators appropriate domains will be nonlinear subspaces of (5.2), and they will be investigated in a follow-up work.

Variational principle for the Vlasov–Poisson equations

In the full Vlasov–Maxwell system, the scalar and vector A potentials are independent dynamic variables, and as such have to appear explicitly in the action functional alongside the stochastic processes X, V and P. In order to ensure the correct coupling between the stochastic processes and the electromagnetic field, an expected value was necessary in the definition of the action functional (5.1). This created a difficulty in deriving a variational principle, as pointed out in the remark following theorem 5.2. This difficulty can be circumvented for the Vlasov–Poisson equations because in this case the electrostatic potential is uniquely determined by the stochastic process X, as will be demonstrated below.

The collisional Vlasov–Poisson equations

The collisional Vlasov–Poisson equations where and and the charge density ρ is given by (1.3), are an approximation of the Vlasov–Maxwell equations in the non-relativistic zero-magnetic field limit. The associated stochastic differential equations take the form and The equations (2.4b), (6.2) and (6.3) form a stochastic reformulation of the Vlasov–Poisson equations. A stochastic particle discretization and the corresponding stochastic variational principle can be derived just like in §§3 and 4, respectively. Also, a variational principle analogous to the Lagrange–d’Alembert principle presented in §5 can be derived in a similar fashion. However, by doing so, one encounters the same difficulty with including the Stratonovich integral. In the case of the Vlasov–Poisson equations a different variational principle can be obtained by observing that the electrostatic potential can be expressed as a functional of the stochastic process X, by solving Poisson’s equation (6.2b). Given the charge density function (2.4b) and specific boundary conditions, the solution of Poisson’s equation can be written using an appropriate Green’s function for the Laplacian. Assuming the spatial domain is unbounded, the standard Green’s function yields From (6.2a) we have the electric field

Let us consider the action functional defined by the formula where the electrostatic potential is given by (6.5). Note that similar to S in (4.4), the functional is itself random, and can be viewed as the action functional of particles represented by the process X which are moving in the electric field generated by particles represented by the process . Similar to the calculations in §4b, the variations of with respect to X, V and P are given by, respectively, where the electric field E is given by (6.6). Note that we are not considering variations with respect to . Let us for convenience define the joint variation of with respect to X, V and P as

In the following theorem, we formulate a variational principle for the system of equations (2.4b), (6.2) and (6.3). Note that is the electric field generated by a distribution of charged particles represented by the process X at time t, and evaluated at the random point in space. Furthermore, the notation means that the variation of is evaluated for the arguments with .

Theorem 6.1 (Stochastic Lagrange–d’Alembert principle for the VP equations).

Let and be stochastic processes, and let be given by (6.5). Assume that and for are functions of their arguments, where f is given by (2.4a). Then , and satisfy the system of stochastic differential equations on the time interval , if and only if they satisfy the following variational principle for arbitrary variations , and , with almost surely, where the action functional is given by (6.8).

Analogous to the proof of theorem 4.2.

It is straightforward to see that equations (6.11), together with (6.5) and (6.6), are equivalent to the system of equations (2.4b), (6.2) and (6.3). The Lagrange–d’Alembert principle (6.12) is unusual in that the variations of the action functional with respect to the argument are omitted. Thanks to such a form, however, the action functional does not require an expected value, and the collisional effects can be correctly included. A similar idea to solve Poisson’s equation and plug the solution into the action functional was presented in [15], where the authors proposed a variational principle for the collisionless Vlasov–Poisson equations. In that approach the energy of the electric field was also included in the variational principle, and the variations were taken with respect to all arguments of the action functional. This approach could be adapted to the stochastic reformulation of the Vlasov–Poisson equations, but the corresponding action functional would have a form similar to (5.1), that is, it would need to contain an expected value, and therefore we would face a similar difficulty as for the Vlasov–Maxwell equations in §5b.

Summary and future work

In this work, we have considered novel stochastic formulations of the collisional Vlasov–Maxwell and Vlasov–Poisson equations, and we have identified new stochastic variational principles underlying these formulations. We have also proposed a stochastic particle method for the Vlasov–Maxwell equations and proved the corresponding stochastic variational principle.

Our work can be extended in several ways. The stochastic variational principle introduced in §4 can be used to construct stochastic variational PIC numerical algorithms for the collisional Vlasov–Maxwell and Vlasov–Poisson equations. Variational integrators are an important class of geometric integrators. This type of numerical scheme is based on discrete variational principles and provides a natural framework for the discretization of Lagrangian systems, including forced, dissipative or constrained ones. These methods have the advantage that they are symplectic when applied to systems without forcing, and in the presence of a symmetry, they satisfy a discrete version of Noether’s theorem. For this reason, they demonstrate superior performance in long-time simulations; see [3,75-84]. Variational integrators were introduced in the context of finite-dimensional mechanical systems, but were later generalized to Lagrangian field theories [2] and applied in many computations, for example in elasticity, electrodynamics, fluid dynamics, or plasma physics; see [10-12,65,72,85-87]. Stochastic variational integrators were first introduced in [35] and further studied in [18,34,37,39,40].

In §5, we have proposed a general action functional for the collisional Vlasov–Maxwell equations. However, we have also determined that in order to prove a relevant variational principle, the domain of this action functional has to be restricted in a way compatible with the collision operator of interest. We have shown that for a class of collision operators with constant diffusion terms, a suitable subdomain is an affine subspace (i.e. a submanifold). A natural continuation of our work would be to investigate submanifolds of (5.2) which are suitable for other collision operators.

Another aspect worth a more detailed investigation is the issue of existence and uniqueness of the solutions of the stochastic reformulations presented in this work, which are non-trivial systems of coupled stochastic and partial differential equations. This question is closely connected to the issue of existence and uniqueness of the solutions of the original collisional Vlasov–Maxwell and the Vlasov–Poisson equations. General results are available in the collisionless case (e.g. [88-90]), but the theory for the collisional equations is less developed (see [29,91-96] and references therein).

Furthermore, our stochastic Lagrange–d’Alembert approach could also be adapted to relativistic plasmas [97], and to variational principles with phase-space Lagrangians appearing in gyrokinetic [13,14,98] and guiding-centre theories [99-102]. In particular, considering stochastic extensions of the variational principles proposed in [99] could offer an alternative stochastic description of anomalous transport in magnetically confined plasmas [103,104].

Finally, as is typical for particle methods in general, the stochastic particle discretization proposed in §3 will require a large number of particles for accurate numerical simulations, which is computationally expensive. Structure-preserving model reduction methods [105,106] have been recently successfully applied to particle discretizations of the collisionless Vlasov equation [107]. It would be of great practical interest to combine our results with model reduction techniques in order to develop new efficient structure-preserving data-driven numerical methods for the collisional Vlasov–Maxwell equations.