Warning: this page is better viewed with Firefox or Safari

*Authors of the code:* S. Colombi & J. Touma

*Publications introducing and exploiting the code* (to be cited when using `VlaPoly`

):

*Obtain* `VlaPoly`

: on request to S. Colombi -> mailto:colombi [snail] iap [period] fr

## Table of Contents

## The waterbag method

In the proper units, Vlasov-Poisson equations are given in one dimension by

where *x* is the position, *v* the velocity, *t* the time, *f(x,v,t)* the phase-space density distribution function, Φ*(x,t)* the gravitational potential and ρ*(x,t)* the projected density.

The waterbag method (DePackh 1962; Robert & Berk 1967; Janin 1971; Cuperman, Harten & Lecar 1971a,b) exploits directly Liouville theorem, namely that the phase-space distribution function is conserved along trajectories of test particles:

This property can be exploited in a powerful way by decomposing the initial distribution on small patches, the *waterbags*, where *f* is approximated by a constant. In one dimension, the numerical implementation is thus conceptually very simple: one just needs to follow the boundaries of the waterbag with a polygon, which is enriched with new vertices when the shape of the waterbag gets more involved (Figure 1).

**Figure 1:** decomposition of phase-space on triangular waterbags and subsequent evolution of the phase-space distribution function. Images: S. Colombi & J. Touma.

## Practical implementation

### Initial conditions

In `VlaPoly`

, the initial phase-space distribution function can be sampled initially either with an adaptive triangulation (Figure 1) or with isocontours of *f(x,v)* (Figure 2). In the isocontour method, which is more optimal, the isocontours are chosen such to bound the mean square difference between the true and the sampled phase-space distribution function weighted by the waterbag thickness. To draw the iscontours, the so-called marching square algorithm is used, inspired from its famous three-dimensional alter-ego (Lorensen & Cline, 1987).

### Poisson equation

The equation of motion of the polygon vertices is the same as for test particles, of which the acceleration is given in one dimension by

where *M*_{right}(*x,t*), *M*_{left}(*x,t*) and *M*_{tot} are respectively the total mass at the right of position *x*, the total mass at the left of position *x* and the total mass.

If one approximates *f* by a constant with value *f*_{k} within a patch, *P*_{k}, *k=*,1...,*N*_{patch}, then

where *s* is a curvilinear coordinate, the last equation making use of Green's theorem. Hence, Poisson equation resolution reduces to a (partial) circulation along the contours of each individual patch that are themselves followed with an adaptive orientated polygon. The calculation of the acceleration corresponds to the most costly part of the code, as many segments of the polygon can contribute to computing the force at a given point when mixing becomes strong.

### Time integration

The vertices of the adaptive polygon are moved with a standard predictor-corrector scheme, with a slowly varying time step, of which the value is mainly determined by the condition

where ρ_{max} is the maximum value of the projected density. This criterion on the time step is derived from harmonic oscillator dynamics considerations. Two other additional constraints are added to make the time step harmonious with refinement of the orientated polygon.

### Refinement

The refinement procedure used in `VlaPoly`

is quite close in spirit to Cuperman et al. (1971a). It consists of a geometric construct using arc of circles passing through sets of three successive points of the polygon and reduces, in the small angle approximation to linearly interpolating local curvature. Refinement criterion relies on local measurement of the surface of the triangle formed by the refinement candidate point and its too neighbors. If the surface of such a triangle exceeds some threshold fixed by the user, then refinement is performed. Additionally, if the distance between two successive points of the polygon exceeds some threshold fixed by the user, a new point is added. There is as well the possibility of unrefining, but this option turns to be suboptimal due to the cumulative noise it introduces in the long run.

### Diagnostics

Diagnostics, in addition to standard conservation tests on energy, volume and mass, include output of snapshots, such as the full polygon structure. This structure can be post-treated to compute various profiles as as well as images of the phase-space distribution function in *(x,v)* or Action-Angle space. To generate these images, we use a modified version of the parity algorithm.

## Illustrations and movies

**Figure 2:** "Landau damping:" evolution of an initially Gaussian phase-space distribution function in *(x,v)* and Action-Angle space to emphasize on the quiescent nature of mixing, as studied in Colombi S., Touma J. (2014). Note: this figure is analogous to Figure 1, but waterbag borders coincide with isocontours of the phase-space distribution function.

- A Movie (same initial conditions, but different simulation): gaussian.avi
**Figure 3:** Mosaic of the evolution in phase-space of systems with various initial conditions. Images: S. Colombi & J. Touma.

## Bibliography

- Cuperman S., Harten A., Lecar M., 1971a, Ap&SS, 13, 411
- Cuperman S., Harten A., Lecar M., 1971b, Ap&SS, 13, 425
- DePackh D. C., 1962, J. Electr. Contrib., 13, 417
- Janin G., 1971, A&A, 11, 188
- Lorensen W.E., Cline H.E., 1987, Comput. Graph., 21, 163
- Roberts K.V., Berk H.L., 1967, Phys. Rev. Lett., 19, 297