## 1 Introduction

The main goal of this letter is to demonstrate the usefulness of symplectic model reduction techniques for numerical simulations of the Vlasov equation using particle methods.

### 1.1 Particle methods for the Vlasov equation

In this work we consider the Vlasov equation

(1.1) |

for the particle density function , where is an external electrostatic field with the potential , and and

are vectors in

with . The standard approach to particle-in-cell (PIC) methods consists of the Ansatz(1.2) |

for the particle density function, where and represent the position and velocity of the -th particle, respectively, and its weight. Substituting (1.2) in (1.1), one obtains a system of ordinary differential equations (ODEs) for and , namely

(1.3) |

It can be easily verified that (1.1) has the form of a Hamiltonian system of equations

(1.4) |

with the Hamiltonian given by

(1.5) |

where and are vectors in .

### 1.2 Geometric integration

The Hamiltonian system (1.1) possesses several characteristic properties. Its flow preserves the Hamiltonian, i.e. the total energy, as well as the canonical symplectic form . The latter property expressed in terms of the standard basis for takes the form of the condition

(1.6) |

where denotes the Jacobi matrix of the flow map , denotes the canonical symplectic matrix defined as

(1.7) |

and is the identity matrix (see, e.g., [4], [5], [6]).

In principle, general purpose numerical schemes for ODEs can be applied to Hamiltonian systems such as (1.1). However, when simulating these systems numerically, it is advisable that the numerical scheme also preserves geometric features such as symplecticy (1.6). Geometric integration of Hamiltonian systems has been thoroughly studied (see [4], [7], [10] and the references therein) and symplectic integrators have been shown to demonstrate superior performance in long-time simulations of such systems, compared to non-symplectic methods. Long-time accuracy and near preservation of the Hamiltonian by symplectic integrators have been rigorously studied using the so-called backward error analysis (see, e.g., [4] and the references therein).

### 1.3 Symplectic model reduction

For the aforementioned reasons it appears desirable to preserve the Hamiltonian structure also in model reduction. In fact, it has been found that preserving the Hamiltonian structure in the construction of the reduced spaces preserves stability [1], which is not guaranteed using non-structure-preserving model reduction techniques. To this end, standard model reduction techniques such as proper orthogonal decomposition have been modified towards the so-called proper symplectic decomposition [9], which does indeed preserve the canonical symplectic structure of many Hamiltonian systems in the reduction procedure. Similarly, greedy algorithms [1] can be used to construct the reduced basis in a Hamiltonian-structure preserving way, and recently also non-orthonormal bases have been considered [3], showing improved efficiency over orthonormal bases.

### 1.4 Outline

The main content of the remainder of this paper is, as follows. In Section 2 we review several model reduction techniques and set the appropriate notation. In Section 3 we present the results of our numerical experiment demonstrating the applicability of model reduction techniques to particle methods for the Vlasov equation. Section 4 contains the summary of our work.

## 2 Model reduction

In this section we briefly review several model reduction techniques and set the notation appropriate for the problem defined in the introduction.

### 2.1 Proper orthogonal decomposition

Proper orthogonal decomposition (POD) is one of the standard model reduction techniques (see [2], [8]). Consider a general ODE

(2.1) |

Equation (1.1) has this form with , , and . When is a very high number, as is typical for particle methods, the system (2.1) becomes very expensive to solve numerically. The main idea of model reduction is to approximate such a high-dimensional dynamical system using a lower-dimensional one that can capture the dominant dynamic properties. Let be an matrix representing empirical data on the system (2.1). For instance, can be a collection of snapshots of a solution of this system,

(2.2) |

at times . These snapshots are calculated for some particular initial conditions or values of parameters that the system (2.1) depends on. A low-rank approximation of

can be done by performing a singular value decomposition (SVD) of

and truncating it after the first largest singular values, that is,(2.3) |

where is the diagonal matrix of the singular values, and are orthogonal matrices, is the diagonal matrix of the first largest singular values, and and are orthogonal matrices constructed by taking the first columns of and , respectively. Let denote a vector in . Substituting in (2.1) yields a reduced ODE for as

(2.4) |

If the singular values of decay sufficiently fast, then one can obtain a good approximation of for such that . Equation (2.4) is then a low-dimensional approximation of (2.1) and can be solved more efficiently. For more details about the POD method we refer the reader to [2]. In the context of particle methods for the Vlasov equation, the reduced model (2.4) allows one to perform numerical computations with a much smaller number of particles. Note, however, that while (1.1) is a Hamiltonian system, there is no guarantee that the reduced model (2.4) will also have that property. In Section 3 we demonstrate that retaining the Hamiltonian structure in the reduced model greatly improves the quality of the numerical solution.

### 2.2 Proper symplectic decomposition

Note that the Hamiltonian system (1.1) can be equivalently written as

(2.5) |

where and . A model reduction technique that retains the symplectic structure of Hamiltonian systems was introduced in [9]. In analogy to POD, this method is called proper symplectic decomposition (PSD). A matrix is called symplectic if it satisfies the condition

(2.6) |

For a symplectic matrix , we can define its symplectic inverse . It is an inverse in the sense that . Let be a vector in . Substituting in (2.5) yields a reduced equation

(2.7) |

which is a lower-dimensional Hamiltonian system with the Hamiltonian . Given a set of empirical data on a Hamiltonian system, the PSD method constructs a symplectic matrix which best approximates that data in a lower-dimensional subspace. We have tested three PSD algorithms, namely the cotangent lift, complex SVD, and greedy algorithms.

#### 2.2.1 Cotangent lift algorithm

This algorithm constructs a symplectic matrix which has the special block diagonal structure

(2.8) |

where is an matrix with orthogonal columns, i.e. . Suppose snapshots of a solution are given as an matrix of the form

(2.9) |

#### 2.2.2 Complex SVD algorithm

By allowing a broader class of symplectic matrices we may get a more optimal approximation. The complex SVD algorithm constructs a symplectic matrix of the form

(2.10) |

where and are matrices satisfying the conditions

(2.11) |

Suppose snapshots of a solution are given as an complex matrix of the form

(2.12) |

where denotes the imaginary unit. The complex SVD of is truncated after the first largest singular values, that is,

(2.13) |

The matrices and are then chosen as the real and imaginary parts of , respectively, that is, . More details can be found in [2].

#### 2.2.3 Greedy algorithm

A greedy approach to the construction of an optimal symplectic matrix has been proposed in [1]. In this algorithm, the matrix is derived column by column via error minimization and symplectic Gram-Schmidt orthogonalization. In contrast to the cotangent lift and complex SVD methods, the greedy algorithm does not assume that the matrix has any specific structure. This typically results in a more accurate reduced system. For more details regarding this algorithm the reader is referred to [1].

## 3 Numerical experiment

In this section we present the results of a simple numerical experiment demonstrating the applicability of model reduction techniques to particle methods for the Vlasov equation.

### 3.1 Initial and boundary conditions

We consider the Vlasov equation (1.1) on a periodic one-dimensional () spatial domain with the initial condition

(3.1) |

where the parameters are set as follows:

(3.2) |

This is a “bump-on-tail” distribution in velocity with a periodic spatial perturbation. The initial conditions for the particle positions and velocities in (1.1

) are generated as random variables drawn from the probability distribution (

3.1) using rejection sampling.### 3.2 Empirical data

We consider a spatially periodic external electric field in (1.1), namely

(3.3) |

where the real parameter is the amplitude. Suppose we have the following computational problem: we would like to scan the domain of , that is, compute the numerical solution of (1.1) for a large number of values of . Given that in practical applications the system (1.1) is very high-dimensional, this task is computationally very intensive. Model reduction can alleviate this substantial computational cost. One can carry out full-scale computations only for a selected number of values of . These data can then be used to identify reduced models, as described in Section 2. The lower-dimensional equations (2.4) or (2.7) can then be solved more efficiently for other values of , thus reducing the overall computational cost. For our simple experiment, we calculated full-scale solutions for the following six values of the parameter :

(3.4) |

The computations were carried out with particles. It should be noted that the same initial values for the positions and velocities of the particles were used in each simulation. The solutions were calculated on the time interval with the time step using the symplectic Störmer-Verlet method. Then, following the description of each algorithm in Section 2, reduced models were derived. The decay of the singular values for the POD, PSD cotangent lift and complex SVD methods is depicted in Figure 3.1.

### 3.3 Reduced model simulations

To test the accuracy of the considered model reduction methods, we compared the results of reduced model simulations to a full-scale reference solution. The reference solution for was calculated on the time interval in the same way as the empirical data in Section 3.2. The POD model (2.4) was solved on the same time interval using the classical 4-th order Runge-Kutta method with the time step for and (thus reducing the dimensionality of the problem from to 10 and 20, respectively). Similarly, the PSD model (2.7) was solved using the symplectic Störmer-Verlet method, where for the cotangent lift and complex SVD algorithms the first and singular value were used, and for the greedy algorithm and basis vectors were calculated, in all cases reducing the dimensionality to 10 and 20, respectively. The relative error of the solutions of the reduced models measured with respect to the reference solution is depicted in Figure 3.2, Figure 3.3, Figure 3.4, and Figure 3.5. We can see that the error for the PSD methods stays bounded on the interval , that is, where the empirical data was available, and starts growing afterwards. The POD solution becomes unstable immediately, which shows how important it is to retain the Hamiltonian structure of the reduced equations.

The total energy (1.5) of the particles for each of the algorithms is depicted in Figure 3.6. We can see that the PSD algorithms preserve the total energy very well. On the other hand, the performance of the POD method is unsatisfactory and it does not appear to improve when more singular values are used. This is a consequence of the fact that the reduced equation (2.4) is not Hamiltonian.

## 4 Summary and future work

We have compared several model reduction techniques and demonstrated their usefulness for particle-based simulations of the Vlasov equation. We have pointed out the importance of retaining the Hamiltonian structure of the equations governing the evolution of particles. Our work can be extended in several directions. First, model reduction methods can be applied to the Vlasov equation coupled to a self-consistent electric field satisfying the Poisson equation, or an electromagnetic field satisfying the Maxwell equations. It would also be interesting to consider collisional Vlasov equations.

## Acknowledgments

We would like to thank Christopher Albert, Eric Sonnendrücker and Udo von Toussaint for useful comments and references. The study is a contribution to the Reduced Complexity Models grant number ZT-I-0010 funded by the Helmholtz Association of German Research Centers.

## References

- [1] B. Afkham and J. Hesthaven. Structure preserving model reduction of parametric Hamiltonian systems. SIAM Journal on Scientific Computing, 39(6):A2616–A2644, 2017.
- [2] A. C. Antoulas, D. C. Sorensen, and S. Gugercin. A survey of model reduction methods for large-scale systems. Contemporary Mathematics, 280:193–219, 2001.
- [3] P. Buchfink, A. Bhatt, and B. Haasdonk. Symplectic model order reduction with non-orthonormal bases. Mathematical and Computational Applications, 24, 2019.
- [4] E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Springer Series in Computational Mathematics. Springer, New York, 2002.
- [5] D. D. Holm, T. Schmah, and C. Stoica. Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions. Oxford Texts in Applied and Engineering Mathematics. Oxford University Press, Oxford, 2009.
- [6] J. Marsden and T. Ratiu. Introduction to Mechanics and Symmetry, volume 17 of Texts in Applied Mathematics. Springer Verlag, 1994.
- [7] R. I. McLachlan and G. R. W. Quispel. Geometric integrators for ODEs. Journal of Physics A: Mathematical and General, 39(19):5251–5285, 2006.
- [8] B. Moore. Principal component analysis in linear systems: Controllability, observability, and model reduction. IEEE Transactions on Automatic Control, 26(1):17–32, 1981.
- [9] L. Peng and K. Mohseni. Symplectic model reduction of Hamiltonian systems. SIAM Journal on Scientific Computing, 38(1):A1–A27, 2016.
- [10] J. M. Sanz-Serna. Symplectic integrators for Hamiltonian problems: an overview. Acta Numerica, 1:243–286, 1992.