Warning: this page is better viewed with Firefox or Safari

*Author of the code:* S. Colombi

*Publication introducing the code* (to be cited when using `Vlamet`

):

*Obtain* `Vlamet`

(current implementation only in 2D phase-space): on request to S. Colombi -> mailto:colombi [snail] iap [period] fr

## Table of Contents

## The semi-Lagrangian approach and the splitting algorithm

We want to follow the evolution of a smooth self-gravitating collisionless fluid following Vlasov-Poisson equations,

where *f(***r**,*u*,*t)* represents the phase-space density at position *r*, velocity *u* and time *t*, Φ is the gravitational potential and *G* is the gravitational constant. We consider the warm case, i.e. the case where the phase-space distribution function is a smooth function of the vector

and can thus be represented on an Eulerian grid.

Most of the numerical methods trying to resolve directly the dynamics of function *f(***P***,t)* come from plasma physics. They are mainly of *semi-Lagrangian* nature and exploit directly Liouville theorem, namely that *f* is conserved along trajectories. In practice, a test particle is associated to each grid site where the value *f*_{T} of the phase-space distribution function has to be evaluated. This test particle is followed backwards in time to find its position at previous time step. Then, function *f* from previous time step is interpolated at the root of the characteristic to obtain the searched value of *f*_{T}, following Liouville theorem.

This approach was pioneered by Cheng & Knorr (1976) in a split fashion along velocity and positions axes. For instance, in the algorithm `VlaSolve`

, we propose an implementation of the splitting algorithm in spherical symmetry using the classical implementation based on third order splines to perform interpolation of the phase-space distribution function.

## The metric method

The main problem common to all implementations of the splitting algorithm is the unavoidable diffusion due to repeated resampling of the phase-space distribution function. Indeed, because of the finite grid resolution, repeated interpolation, as sophisticated it can be, induces information loss and augmentation of entropy.

The way we proceed in `Vlamet`

to improve upon this issue consists in employing metric elements to follow closely the flow and its deformation during several time steps, hence reducing considerably diffusion. Namely, these metric elements, of which the position *P*_{m}*(***Q**_{m},t) follows the Lagrangian equations of motion (*Q*_{m} being the initial, Lagrangian position at *t=0*), transport as well numerically the deformation tensor *T*_{m}

and the Hessian of the displacement *H*_{m}

This is symbolically illustrated in Fig. 1 by the thick lines which intersect at metric element positions, such as in the blue, the red, the orange and the purple points.

To reconstruct the phase-space distribution function *f(***P***,t)* at a point *P* of the grey Eulerian grid of Fig. 1, it is very easy to find the initial (Lagrangian) position *Q**(***P***)* of an arbitrary test particle at (Eulerian) position *P*, by using a second-order approximation of the local metric:

near the closest metric element to the point *P* (the blue point on Fig. 1).

**Figure 2:** Phase-space of a simulation with an initial Gaussian phase-space distribution function (left panel) and the region of influence of each metric element obtained from a percolation algorithm in Lagrangian space (right panel), as displayed with a random colour. From Colombi & Alard (2017).
**Figure 3:** Projected density in a simulation with Gaussian initial conditions. Simulations performed with the metric approach (I, X, XI) are compared to the results obtained with the splitting algorithm (XII, XIV, XVI), as well as an entropy conserving waterbag run of Colombi & Touma (2014) performed with `VlaPoly`

. The waterbag run is so accurate that he can be considered as representing the "exact" solution. As indicated on the panel, Δ is the size of the step of the grid used to represent the phase-space distribution function, while Δ_{m} indicates the inter-element of metric distance. One can see that in the particular case considered here, the metric approach provides comparable results to splitting algorithm runs with twice the spatial resolution. From Colombi & Alard (2017).
To reach best accuracy, the distance between each metric element and the test particle is measured in Lagrangian space (left panel of Fig. 1), which requires a special percolation algorithm of which the result is illustrated by Fig. 2. Then, once one has access to the initial position *Q* of a test particle, one simply interpolates the initial phase-space distribution function to reconstruct *f(***P***,t)* using Liouville theorem. This is performed quickly using a cheap second order interpolation. At some point it is necessary to create new initial conditions with new isotropic metric elements and repeat the procedure until next resampling. During this step, a third order spline interpolation is used, and *f(***P***,t)* reconstructed this way is taken as a new initial condition, exactly as in the traditional splitting method. Additionally, interpolation between nearby Lagrangian candidates (the blue, red, orange and purple points on left panel of Fig. 1) is performed to preserve smoothness of the procedure.

The metric algorithm is second order in space and time and mathematical convergence of a simplified version of it has been studied in Campos Pinto & Charles (2016).

While our actual implementation is performed in the simplest case, the 2D dimensional phase-space, we discuss in Colombi & Alard (2017) how it can be rather straightforwardly generalised to 4 and 6 dimensional phase-space in a parallel fashion. Additionally, it is obviously generalisable to plasmas.

## Performances

Because the metric method is not of splitting nature, i.e. does not allow one to reduce the phase-space reconstruction to simple one dimensional interpolations on the phase-space grid, it is more expensive than the traditional splitting algorithm at the moments of remap. This is expected to be particularly true in 4D and 6D dimensional phase-space. However, this is compensated by the fact that remappings are seldom performed and that this algorithm is much less diffusive than its splitting counterpart, as illustrated by Fig. 3. In Colombi & Alard (2017) we argue how implementations of `Vlamet`

in 4 and 6D would in fact be highly competitive.

## Bibliography

- Campos Pinto M., Charles F., 2016, preprint HAL-01385676
- Cheng C. Z., Knorr G., 1976, JCoPh 22, 330
- Colombi S., Touma J., 2014, MNRAS 441, 2414