Relativistic Vlasov-Maxwell modelling using finite volumes and adaptive mesh refinement

The dynamics of collisionless plasmas can be modelled by the Vlasov-Maxwell system of equations. An Eulerian approach is needed to accurately describe processes that are governed by high energy tails in the distribution function, but is of limited efficiency for high dimensional problems. The use of an adaptive mesh can reduce the scaling of the computational cost with the dimension of the problem. Here, we present a relativistic Eulerian Vlasov-Maxwell solver with block-structured adaptive mesh refinement in one spatial and one momentum dimension. The discretization of the Vlasov equation is based on a high-order finite volume method. A flux corrected transport algorithm is applied to limit spurious oscillations and ensure the physical character of the distribution function. We demonstrate a speed-up by a factor of 7 × in a typical scenario involving laser pulse interaction with an underdense plasma due to the use of an adaptive mesh.


Introduction
The Vlasov-Maxwell system of equations describes the dynamics of a collisionless plasma whose component species interact through self-consistent electromagnetic fields.It is of critical importance to the fundamental understanding of non-equilibrium processes in collisionless plasmas, as well as many practical applications, for example laser plasma acceleration [1,2], inertial confinement fusion [3][4][5], high harmonic generation [6] and shocks in astrophysical plasma [7].
Numerical approaches to solve this system are primarily divided into Particle-In-Cell (PIC) methods, which approximate the plasma by a finite number of macro-particles, and methods that discretize the distribution function on a grid: so-called Eulerian methods [8].As PIC methods do not require a grid in momentum space, they are efficient at handling the large range of scales associated with relativistic laser-plasma interaction.They are therefore suitable for modelling high dimensional problems [9].However, the approximation of the distribution function by a finite number of particles introduces statistical noise, making it difficult to resolve fine velocity space structures and high energy tails within the distribution function.
The ability to resolve fine structures related to low density tails in the distribution function is of critical importance to topics in laser plasma acceleration, e.g., modelling of collisionless shock acceleration (CSA).In laser plasma induced CSA, a small fraction of the ion population is reflected by an electrostatic potential barrier set up by laser plasma interaction.This low density tail of the distribution function has been suggested to play an important role for the dynamics [10], indicating that high resolution is needed.
Eulerian methods, which discretize the distribution function on a grid have very low levels of numerical noise.They are therefore appropriate for the detailed study of processes where a small number of high energy particles play a significant role.The most widely used method for solution of the Vlasov equation is time-splitting, first suggested by Knorr and Cheng [11].The method involves splitting the Vlasov equation into lower dimensional advection equations that are alternately advanced.Second order accurate time splitting methods have been used to solve the Vlasov-Poisson and Vlasov-Maxwell system of equations in Refs.[11][12][13][14][15][16][17].In these references the distribution function is represented on a fixed uniform grid, which leads to large computational costs, especially for multidimensional simulations.
To increase the computational efficiency and be able to treat problems with a wide range of scales, the Vlasov-Maxwell system can be represented on an adaptive mesh in combination with higher order methods.For this approach, high-order finite volume discretizations can be used to solve the Vlasov-Maxwell equation system, see for example Ref. [18].The adaptive mesh then evolves as the characteristics of the distribution function develop, which allows higher resolution to be applied to those parts of phase space that exhibit complicated behaviour.At the same time, the distribution function remains well-resolved in regions with a coarse mesh through the use of high-order numerical schemes.This means that the use of adaptive meshing limits the computational effort to regions with small scales but still maintains a high degree of accuracy in the full domain.
Adaptive Mesh Refinement (AMR) is not a new concept, and has been used extensively to minimise computational overhead in a variety of systems.Some examples focusing on laser plasma research include: the PIC code Warp [19], which implements AMR for its Poisson solver; M. Dorr et al. [20] investigates speckles and filamentation in inertial confinement fusion by solving the Poisson-Euler equations coupled to Maxwell equations; Ref. [21] which solves the Boltzmann equations using a hybrid octree-AMR approach.AMR has also been used in Ref. [22] to simulate the Vlasov-Poisson system, i.e. the electrostatic and non-relativistic limit of the Vlasov-Maxwell equations.Finally, N. Besse et al. [9] use a wavelet based adaptive grid to solve the relativistic Vlasov-Maxwell system.
In this paper we present the open source blockstructured Eulerian Vlasov-Maxwell solver veritas (Vlasov EuleRIan Tool for Acceleration Studies) [23].The solver is based on a high-order finite volume method, implementing the flux corrected transport algorithm to limit spurious oscillations in the distribution functions, in the presence of steep gradients.veritas offers the capability to study realistic laser-plasma problems in two dimensions (1D1P), with a complete electrodynamic framework (i.e.Vlasov-Maxwell) at relativistic speeds.To our knowledge, veritas is the first complete, relativistic Vlasov-Maxwell solver using AMR which shows a significant performance increase.This advancement moves continuum solvers towards the category of capable simulation tools alongside their PIC counterparts.
The paper is organized as follows.Section II describes the Vlasov-Maxwell equations.Section III presents the numerical scheme for veritas, including a description of the mesh structure, information flow, regridding procedure, finite volume discretization of the Vlasov equation and the discretization of Maxwells equations.Section IV describes benchmarking with comparison to results from analytical theory and PIC simulations, and demonstrates the improved performance of the adaptive mesh approach.Conclusions are summarized in Section V.

The Vlasov-Maxwell system
The Vlasov-Maxwell equation system describes the time evolution of the electron and ion distribution functions, which interact self-consistently with the electromagnetic fields in a collisionless plasma.For the case of a plasma with spatial variation in one direction, the Vlasov-Maxwell equations can be reduced to a two dimensional 1D1P problem: where f s is the distribution function of a species (e.g.electrons or ions), E and B are the electric and magnetic fields, x is a spatial coordinate, p x is a momentum coordinate in the direction of x, q is the charge, m denotes the rest mass of the charged particles (electrons or ions) and γ = p 2 /m 2 c 2 + 1 is the relativistic factor.The singleparticle Hamiltonian yields conservation relations for the transverse canonical momentum (orthogonal to the direction of variation of the plasma): Π ⊥ = qA ⊥ + p ⊥ = 0 [24,25].The conservation of Π ⊥ stems from the fact that the perpendicular coordinates y and z do not enter the Hamiltonian.Here, c is the speed of light and φ and A are the electrostatic and vector potentials, respectively.For a one-dimensional system, Maxwell's equations take the form The distribution function is represented on a block structured mesh, which is adapted during time evolution to ensure that regions of phase space with more complex dynamics are associated with a finer resolution in the mesh [26].We use a finite volume scheme which is fourthorder accurate, making it possible to use a very coarse representation of the distribution function in regions with less complex dynamics, without introducing numerical instabilities.
The electromagnetic fields are defined on a mesh associated with the finest spatial resolution of the distribution function, using a fourth-order discretization of Maxwell's equations and a staggered grid for the electric and magnetic fields.As the fields are one dimensional, the use of a mesh with the finest spatial resolution as opposed to using adaption to adjust the resolution comes at only a moderate cost.
Distribution functions, for one or multiple plasma charge species, are time advanced together with the electromagnetic fields, yielding an overall self-consistent solution to the Vlasov-Maxwell system.In the following, we describe the mesh structure, the representation of the solution, the adaption and information flow from one mesh to another, as well as the discretizations used in veritas.

Mesh structure
To resolve the different scales of the distribution function f (x, p x ), the domain D = [x min , x max ] × [p x,min , p x,max ] is discretized using a block structured mesh.The mesh consists of a number of levels L i of different resolution, where i = 0, . . ., n max , ranging from coarsest to finest.Each level is a union of rectangular patches, i.e.L i = ∪ j R i,j , where R i,j denotes a rectangular patch on level i.Each rectangular patch, which is disjoint from all other rectangular patches on the same level, consists of a number of rectangular cells of side lengths ∆x i = ∆x 0 /r i and ∆p x,i = ∆p x,0 /r i , where r is the refinement ratio.
The distribution function is described by cell-averaged values f m,n k,l : where m, n are indices for the level of refinement along the spatial and momentum dimension, respectively, and are the bounding dimensions of the cell.Throughout this paper we indicate a cell-averaged value with a tilde (˜).Concerning the representation of the distribution function on the block structured mesh, the indices m and n always take the same value.The use of cells in the mesh with different m and n would correspond to independently adapting in the x and p x direction.Although in principle possible, this would complicate the mesh structure and is currently forbidden.On the other hand, in the calculation of the current and charge densities, the distribution function is interpolated to the finest spatial mesh; which may result in values in cells with different indices for m and n.
The coarsest level L 0 consists of a single rectangle and finer levels are nested in coarser levels L i+1 ⊂ L i .Figure 1 shows a rectangular patch on level L i+1 which is partially overlapped by a rectangular patch at level L i .The finer rectangular patch is shaded blue with a blue dot at the center of each cell and the coarser cells are indicated by black dots at the cell centers.Blue dots outside the blue shaded region represent ghost cells which are used to time advance the interior cells.Ghost cells are interpolated from the cell-centered values in the overlaid coarser level, as will be described in Section 3.3, or read from an adjacent rectangle on the same level (i.e.interior ghost cells).Values of the distribution function on a coarser level which is overlaid by a finer level, e.g. the black dotted cells inside the blue shaded region, are defined by the average of the values on the finer level f i,i k,l = 1/r the value for the ghost cells which are interior to the level.3. From the coarsest to the finest level, for each pair of levels L i , L i+1 , interpolate ghost cells for rectangles in L i+1 which are not interior, using values from L i .
The time stepping procedure itself is outlined in Section 3.7.

Adaptive mesh refinement
Figure 2 shows a plasma which has interacted with a laser pulse.The distribution function is represented on five levels which are divided into rectangles, with high resolution in regions with higher densities and more complex dynamics.To track these regions, the mesh structure is updated in the following way: 1.For each level L i , mark each cell that belongs to a rectangle in the level and has an error indicator exceeding some threshold δ thres , discussed in Section 3.8.The marked points M i are used to generate a new level and rectangle structure L i where i = 0, ..., n max .2. Let n max − 1 be the largest integer such that M n max −1 is non-empty and create a new level L ñmax such that Furthermore, upon adapting the mesh, the representation of the solution on the old mesh structure must be transformed into a representation of the solution on the new one.Information contained in a cell in the new level L i may originate from L i or a subset of a cell in L j for some j < i.In the former case, the cell-averaged value in L i is taken directly from L i , whereas in the latter case an interpolation is performed.This operation is performed for all L i before interpolating those values that cannot be copied from L i .As the coarsest level is the same in both the old and new mesh, all values in L 0 are determined from values in L 0 .Values in cells of L i that are not a subset of L i can then be interpolated from L i−1 for i = 1, ..., n max .
The covering of M i ≡ L i+1 ∪ M i−1 with rectangles closely follows the AMR implementation pioneered by Berger et al. [26,27], with specific extensions outlined in their follow-up papers [28,29].To keep dependencies to a minimum, and integration tight & efficient, no external libraries are called here and the AMR infrastructure outlined below has been implemented directly into veritas.For each level, from the finest to coarsest, we identify a minimum bounding box R i,bb , such that all marked cells on level i are extant within the rectangle's boundaries.R i,bb is then split into smaller rectangles R i,j through a recursive process until a minimum efficiency, defined as the proportion of marked cells in each rectangle, is reached: Here, N [•] denotes the number of cells within a given set, and ε min is a user-set parameter for the minimum efficiency.The specific value of the parameter will affect the overall runtime of the program.A low value for ε min leads to a larger mesh with more degrees of freedom, with a corresponding increase in computational work.On the other hand, a high value leads to a more complicated mesh structure, with larger overheads, which in addition must be updated more frequently to accurately represent the solution.Values for ε min between 60 % and 80 %, are found to give reasonable performance, although optimal values are problem dependent.
To split a rectangle, we introduce signatures Σ x,k and Σ px,l , which are functions of the discrete x and p x coordinates, respectively.For a given rectangle, Σ x,k is defined as the number of cells in the intersection of the rectangle and marked cells M i that have the spatial index k.Σ px,l is defined similarly, but with the spatial and momentum coordinates exchanged.A rectangle can be split into two smaller rectangles at an index, either k or l, at which the corresponding signature is zero.If neither signatures contain a zero, the signature's derivatives ∆ x and ∆ px are used, and rectangle edges are identified where zero crossings occur.In the case of more than one zero crossing, the crossing with largest rate of change in ∆ is chosen for the partition; if two crossings have the same magnitude, the one closest to the rectangle center is chosen to prevent thin rectangles which reduce efficiency.If none of the above partitioning criteria are met and the efficiency ratio still has not satisfied the ε min value, the rectangle is bisected along its longest dimension.

Coarse to fine interpolation
High-order interpolation of the distribution function from a coarser level to a finer level is performed (1) to calculate cell-averaged values in ghost-cells, (2) to interpolate charge and current densities to the finest level of the grid and (3) to interpolate data to a refined cell after performing mesh adaption.For high-order discretizations, which are needed (for example) to treat the disparate scales of the discretization in an adaptive solver, low order slope limited interpolation is not suitable.Instead methods such as filtered high order interpolations [30], Weighted Essentially Non-Oscillatory (WENO) techniques [31][32][33] and least squares methods [34] can be used.
In this work, we use a simple fourth-order conservative least square interpolation method.Coarse to fine interpolation of the distribution function, i.e. interpolation from f i,i k,l to f i+1,i+1 rk+m,rl+n , where m, n = 0, . . ., r − 1, is performed in two steps, by first interpolating f i,i k,l to f i+1,i rk+m,l and then to f i+1,i+1 rk+m,rl+n .Hence, it is sufficient to describe only the interpolation step from f i,i k,l to f i+1,i rk+m,l , because the interpolation in the other coordinate is analogous.This is performed by first introducing which is interpolated using a third degree polynomial ), ensuring particle number conservation.The values for f i+1,i rk+m,l are calculated from ), where x i+1 rk+m± 1 2 are the end-points of the cells in the refined level.

Discretization of the Vlasov equation
We introduce a finite volume method which decomposes the discretized advection operator into fluxes across cell boundaries.On each level i in the mesh, the Vlasov equation is averaged over cells in the mesh, resulting in a set of ordinary differential equations (ODEs): where f i,i k,l is the cell-averaged value of the distribution function, and denote fluxes.The fluxes are defined as: and where The ODEs ( 9) are exact and the numerical approximations enter in the calculation of the flux terms and . Despite the gain from the fast convergence of high-order methods in the calculation of the flux terms, high-order methods suffer from spurious oscillations in regions with under-resolved gradients that could trigger numerical instabilities.To avoid this, we use a Flux Corrected Transport (FCT) algorithm in the calculation of the advective terms in the Vlasov equation as suggested in Ref. [18].The implementation of FCT, which is described in 3.7, mixes low and high-order fluxes in such a way as to maximize the high-order flux without introducing unphysical properties in the distribution function.Whilst investigations of modified and refined FCT methods to limit spurious oscillations associated with high-order methods are underway [35], we use the original implementation in Ref. [36].The reduction in order of accuracy for the limited solution in regions with under-resolved gradients is compensated by the use of a finer grid in exactly these regions.
We evaluate each flux in two different ways, using a low and a high-order method, respectively.The low order fluxes, denoted , are evaluated using first order upwinding, i.e. as the product of the face-averaged force terms and the cell-averaged value of the distribution function in the upwind cell.In section 3.7, the advection operator is calculated by blending these with high-order fluxes , using the FCT algorithm, to obtain a stable scheme which does not create unphysical extrema in the distribution function.
To obtain a second order accurate calculation of , the cell-averaged and faceaveraged quantities can be approximated by their cell and face centered values.However, for a finite volume scheme beyond second order accuracy it is necessary to distinguish cell-centered and averaged values as well as using accurate quadrature rules for the face-integrals.Here, we follow Ref. [37] and add corrections to the midpoint approximation of the face-integrals using the transverse derivatives, yielding: and respectively.This is a fourth-order accurate expression for the flux-integrals.Here, the face-averaged values of F x and F px are defined as and the distribution functions ) dx, (17) which are used to evaluate the fluxes.
The Hamiltonian for a particle is Using Hamilton's equations ∂H/∂p x = F x and −∂H/∂x = F px , the face-averaged force terms take the forms and This indicates that the vector potential must be evaluated on the spatial faces.For the electric field, we can directly evaluate its cell-averaged value through the electrostatic potential for which we solve.
To evaluate f i k+ 1   2 ,l (and similarly for ) an upwind biased WENO-type reconstruction is used [18,31].The WENO scheme makes use of a four cell stencil involving the nearest and next nearest neighbours in the normal direction of the (k + 1  2 )-face.These four cells are divided into two sub-stencils, defining two third order interpolants and A weighted interpolant is then obtained by setting where 2 ,l .Following the WENO procedure, β L , β R are calculated based on estimates for the smoothness of the interpolant and approach β L = β R = 1/2 in the limit of a smooth distribution function.However, in contrast to the conventional WENO algorithm [31], we follow Ref. [18], so that the largest weight is assigned to the value of p L or p R that is associated with the upwind stencil.
Up to this point, we have not considered the effect of coarse-fine interfaces in the calculation of fluxes.Coarse to fine interfaces are illustrated in Figure 1.At the cell faces which constitute a coarse-fine interface, i.e. the exterior boundary of a fine level L i+1 , fluxes on the fine level and coarser level L i must be defined consistently in order to obtain conservation of particle number.Conservation up to machine error is critical due to the feedback of the charge density through the longitudinal electric field.For a spatial face (k + 1/2, l) in level L i , the flux is defined using the fluxes calculated on the r constituent finer cell faces: This enforces that the number of particles flowing from the fine cells is the same as the number of particles that enter the coarse cell.Faces associated with the momentum coordinate are treated similarly.

Evaluation of charge and current densities
To evaluate a general moment a reduction operation is performed over the different levels and rectangles, yielding a sum of the form: Here, ] or is otherwise zero, and I(x, l, i) approximates using solution quantities at the level L i .
To evaluate the cell-averaged charge density, values of the distribution function, represented on different levels i, are interpolated using least-square interpolation to the finest spatial level n max .The interpolated values are denoted f nmax,i k,l and the charge density is then calculated by the reduction: The only source of error introduced in this relation is due to the interpolation to the finest grid with respect to the spatial coordinate.Accurate calculation of the current density is more challenging as the integrand in Eq. ( 3) has a more complicated dependence on x and p x (through γ(p x , A ⊥ (x)), f (x, p x ) and A ⊥ (x)), but is of critical importance for the interaction between the plasma and the electromagnetic field.Using the interpolated values f nmax,i k,l , the spatial cell averaging is made to second-order accuracy by commuting the cell averaging and evaluation of the integrand with respect to the spatial coordinate: As for the charge density, Lk is evaluated using Eq. ( 26), but in contrast to Eq. ( 28), the integrand I( Ã⊥,k , l, i) needs further consideration and use of a quadrature formula based on cell averaged quantities.The situation is similar to the flux calculation in 3.4 and the fourth order scheme: is used.In this expression, the momentum average values for 1/γ are calculated analytically:

Discretization of Maxwell's equations
To solve the Vlasov-Maxwell system self-consistently, we perform a spatial discretization of the equations for the transverse fields, resulting in a set of ODEs which are simultaneously time stepped with the kinetic equation.
We discretize E y , E z , A y and A z using the cell-centred values Ẽy,k , Ẽz,k , Ãy,k and Ãz,k , i.e. averages over the same spatial cells as for the kinetic equation.To avoid oddeven decoupling, we represent B y and B z as cell-centred averages By,k+ 1 2 and Bz,k+ 1 2 , at a staggered grid.The fields are discretized on the finest level for the distribution function.This is motivated by the fact that Maxwell's equations only depend on the spatial coordinate making the computation cheap compared to the computation for the Vlasov equations, in spite of the use of a fine mesh.
The equations for the transverse fields are cell averaged and the spatial derivatives are discretized to fourth order, yielding ∂ By,k+ and Here, Jk denotes the cell-averaged current density.In simulations of laser matter interaction, a laser pulse is implemented as a Dirichlet boundary condition at the left hand side of the simulation box.
Regarding the use of the transverse fields to evaluate the coefficients in the Vlasov equation, we note that Eq. (19) and Eq. ( 20) depend on the point values of the transverse vector potential at x nmax k+ 1 2 .To calculate the vector potential on the spatial faces from its cell-averaged values, we use a fourth order WENO interpolation, similar to the one that was used for face interpolation of the distribution function, but in this case without upwind biasing of the smoothness indicators.Including a non-linear scheme here is primarily done for the sake of safety: to squelch numerical sources of noise.Being a one dimensional interpolation, the computational overhead is minimal.
For the equation with the longitudinal electric field E x , we introduce φ, such that E x = −∂φ/∂x.The potential φ satisfies the Poisson equation ∆φ = −ρ/ 0 .A fourth order discretization is given by Ref. [22]: where φk is the cell-averaged potential and ρk is the cellaveraged charge density.
For simulations of laser matter interaction, where no charge leaves the simulation box, we take the electric field at both boundaries of the simulation box to be zero.The boundary conditions are implemented by splitting the potential in two parts φ = φ 1 + φ 2 where φ 1 satisfies Eq. ( 38) with periodic boundary conditions and hence makes it possible to avoid modification of the discretization at the boundary points.Furthermore, we take φ 2 = −E 0 x, which satisfies the homogeneous Poisson equation and choose E 0 such that the homogeneous Neumann condition is fulfilled at the right boundary.
The cell-averaged electric field, which is used to evaluate F px , is then obtained to fourth order from:

Time advancement
For time advancement of the ODEs in the discretization of the Vlasov-Maxwell system we use the Runge-Kutta method ARK4(3)6L [2]SA, which is a six stage, fourth-order accurate method [38].It is a suitable choice for future extensions involving a Fokker-Planck diffusion operator.In such extensions the advection operator would be treated explicitly and the diffusion term implicitly [37].
The combined set of ODEs for the transverse electromagnetic fields and the distribution function can be regarded as a system of equations of the form where L stands for the respective right hand sides of Eqs. ( 9), ( 33)- (37).Denoting the solution at time step t n by F (t n ), a sequence of solutions, F i (t i n ) = F (t n ) + ∆t n L i , at intermediate times t i n = c i ∆t + t n are constructed, for i = 1, . . ., 6. Here, L i is a linear combination of L(F ) at the previous intermediate time steps: Once L F i (t i n ) has been calculated for i = 1, . . ., 6, the solution at time step t n+1 = t n + ∆t is calculated from where the values of a ij , b j and c i can be found in Ref. [38].
The part of L(F ) that corresponds to Maxwell's equations is evaluated using the discretization in subsection 3.6.On the other hand, for the Vlasov part, the FCT algorithm is applied in order to limit anti-diffusive fluxes, which may cause unphysical extrema and is the reason behind introducing the low-order fluxes in subsection 3.4.Although the FCT algorithm is found to be necessary to obtain a stable scheme, its use comes at the expense of introducing a source of hyper-diffusivity in regions with under-resolved gradients.To mitigate this, the criterion for refinement in subsection 3.8 involves terms proportional to derivatives of the distribution function.
When applying the Runge-Kutta method, the flux that is used to calculate the solution at an intermediate or final state in the time stepping algorithm is a linear combination of the fluxes calculated from the intermediate solutions at earlier stages.The generic structure of an advancement of the distribution function takes the form where f 0 k,l , f 1 k,l is the solution before and after time advancement, respectively, and F x , F px denote linear combinations of flux terms calculated from intermediate solutions.Here, we have simplified the notation by ignoring the level in the mesh, on which the quantities are defined.The combined flux terms F x , F px can be evaluated using either the low or high order fluxes, denoted F L x , F L px and F H x , F H px , respectively.For the high order flux, the weights of the fluxes at intermediate steps are those according to the Runge-Kutta method.However, to ensure that the low order solution is positive, the upwind flux at time t n is used.
Following Ref. [36], we first calculate an approximation to f 1 k,l using the low order fluxes This is followed by defining anti-diffusive fluxes and Based on these, the total amount of anti-diffusive fluxes into (P + k,l ) and out (P − k,l ) of a cell is defined by and Furthermore, for cells that share faces with the cell with indices k and l, we define f max k,l and f min k,l to be the maximum and minimum values, respectively, of both f 1,L k,l and f 0 k,l .The maxima and minima are used to define , which describe the maximum allowed increase/decrease that is compatible with not creating unphysical extrema.Based on this, a maximum fraction of the anti-diffusive flux that can be allowed to enter/exit a cell is given by Depending on the sign of the anti-diffusive flux across the different faces, the fraction of the anti-diffusive flux that can be admitted is given by and These fractions are finally used to add the maximum allowed contribution from the higher order fluxes For an explicit time-stepping scheme to be stable, the time step must fulfill a Courant-Friedrichs-Lewy (CFL) condition.Here, the time step is determined by enforcing a user provided maximum CFL number such that the worst case coefficients (which scale with a 0 ) in the simulation do not exceed this value.The available range of CFL numbers is dictated by the stability region of the Runge-Kutta algorithm.

Refinement indicator
To concentrate finer resolution to parts of phase-space where the distribution function has a complicated structure, we choose a refinement indicator which depends on the magnitude of the distribution function as well as its first and second derivatives: Here, the constants w x,1 r −i/2 , w px,1 r −i/2 , w x,2 r −2i , w px,2 r −2i , w f r −i are chosen such that they are all inversely proportional to f max , where f max is the maximum value of the distribution function, which is determined from the initial condition.
In addition to the criteria for mesh refinement, we have implemented the possibility to specify a minimum level of refinement as a function of the phase space coordinates.This makes it possible to guarantee that parts of phasespace which are known to carry important information about the distribution function, and have complex dynamics a priori, are sufficiently resolved.The plasma-vacuum interface is an example of such a region, where most of the laser-plasma interaction occurs.Comparison of analytical solutions for density (dashed lines) and squared normalized vector potential (solid lines) and solutions that are calculated using veritas (black markers), for a circularly polarized laser pulse that impinges on an overdense plasma.The density has been normalized by 2nc, to fit on the same scale as the squared vector potential.The analytical results for different intensities are labelled with the colours red (a0 = 0.25), blue (a0 = 0.50), green (a0 = 0.75) and purple (a0 = 1.00).

Numerical benchmarking
To benchmark veritas, we present results from Vlasov-Maxwell simulations in two cases.In the first case we consider a situation with a circularly polarized (CP) laser pulse impinging on an overdense electron plasma.In this case, it is possible to derive analytical solutions for the electron density by using fluid theory, and these analytical results in turn can be compared with numerical results of veritas.
The second benchmarking example will be a study of laser plasma interaction in three different regimes of laserplasma interaction.Here we will make a comparison with results of the well-established PIC code epoch [39].

Circularly polarized light impinging on a plasma slab -comparison to analytic theory
We study a circularly polarized laser pulse interacting with an overdense plasma, i.e. a plasma with a density n 0 higher than the critical density n c = ω 2 m e 0 /e 2 , where ω is the laser frequency.In particular, we are interested in comparing analytical and numerical solutions in the regime of total reflection for a cold electron plasma (T e m e c 2 ) under the assumption of immobile ions.In this regime a standing wave is formed by the interference of the incoming and reflected pulses, while penetration of the pulse into the plasma is limited to the skin depth.
Using cold fluid theory, analytical solutions can be derived that describe the quasistationary state reached by the plasma and electromagnetic fields [40][41][42].Here we follow the notation of Ref. [43].
In our simulations, we consider a laser pulse incident from vacuum (x < 0) on a semi-infinite plasma slab (x > 0).The initial distribution function of the electrons is where Θ(x) is the Heaviside function.We introduce the normalized vector potential of the incoming laser pulse to be a L (x, t) = |q e |A(x, t)/m e c.For circularly polarized light where ŷ and ẑ denote unit vectors forming an orthonormal basis in the plane transverse to the laser propagation direction.The envelope is where T is the laser period and a 0 is the incident laser field amplitude.The laser pulse is imposed as a boundary condition for the transverse magnetic field on the left side of the simulation box, which is obtained by taking the derivatives B = ∇ × A.
The analytical theory is based on the assumption that the plasma-vacuum interface is pushed up to a point x b where the ponderomotive force −m e c 2 ∂γ/∂x of the laser pulse is balanced by the electrostatic field q e E x due to charge separation.Assuming that an equilibrium has been reached, p x = 0 and the relation between the Lorentz factor and the normalized laser amplitude is γ(x) = 1 + a 2 (x).In the following, we express length in the inverse wave number k −1 = c/ω and density in terms of the critical density.Letting a b denote the value of the vector potential (envelope) at the equilibrium point x b , it can be shown that 2a This relation can be solved numerically for a b and the value x b is then calculated from Using these values for x b and a b , the normalized vector potential in the vacuum and plasma regions become and respectively, where λ s = 1/ √ n 0 − 1 is the skin depth and x 0 is determined by ensuring the continuity of the vector potential at x = x b .The electron density is finally calculated from n e = n 0 + ∂ 2 γ/∂x 2 [43].
Figure 3 shows, that the analytical solutions for the electron density and vector potential given above agree well with solutions obtained using veritas, for a plasma of density n 0 = 2n c and the four different laser intensities a 0 = 0.25, 0.50, 0.75 and 1.00.In the numerical simulations, the semi-infinite cold plasma was represented by a finite plasma slab of length 5λ.The dimensions of the simulation box were [−3, 7]λ × [−8, 8]m e c.The mesh consisted of five levels with refinement ratio r = 2 and, on the coarsest level, n x = 76, n p = 50 points in the spatial and momentum direction respectively.The temperature of the electrons was taken to be T e = 5 • 10 −4 m e c 2 .The finite temperature and additional heating due to the finite rise time of the laser pulse results in a slightly less peaked structure in the electron density at the vacuum plasma interface compared to analytical theory which neglects these effects.However, reducing the temperature significantly would increase the number of points in momentum space that are needed to resolve the p x dependence of the distribution function.

Comparison to PIC simulations in different interaction regimes
We performed simulations in order to study the performance of our code in three different interaction regimes: underdense plasma, the relativistic self-induced transparency (SIT) regime and the hole-boring regime.For relativistic laser pulses with a 0 1 propagating in infinite plasma the classical critical density n c = ω 2 m e 0 /e 2 , introduced in Sect.4.1 has to be modified in order to take into account the dependence of the electron effective mass on the γ-factor.Using the normalizations of Sect.4.1 for CP pulses and taking into account conservation of canonical momentum, this leads to the relativistic critical density n eff c = 1 + a 2 0 /2 n c [44,45].The possibility for propagation of a relativistic pulse in a classically overdense plasma is known as self-induced transparency.However, as we have seen in Sect.4.1, for semi-infinite plasma -a local density peak is formed at the plasma-vacuum interface, leading to a critical density departing from n eff c [41,42].The situation is complicated by kinetic effects [43] and ion motion [46,47] and the threshold for transition from the transparent to the opaque regime has to be determined numerically.In particular, when ion motion effects are taken into account one studies the transition between SIT and the so-called hole-boring regime [48][49][50][51][52] in which the ions are accelerated in the charge-separation-induced electrostatic field and the whole plasma-vacuum interface recedes deeper into the plasma.
We performed simulations in order to study the interaction of plasmas with different densities and a laser pulse with normalized laser intensity a 0 = 2. Figure 4 shows the distribution function for electrons at time t = 15T after the front of the laser pulse has reached the plasma, for the densities n e = 0.3n c (in the underdense plasma regime),

Efficient modelling of laser plasma interaction
The performance of the adaptive solver veritas is evaluated by comparing to the performance of using a uniform grid, i.e. when the adaptive solver is operated with a single level.To identify how the performance depends on the problem, we consider the three cases n e = 0.3n c , 1.5n c and 2.6n c , which were described in subsection 4.2.The case with lowest density has the most complex electron dynamics and will hence be the most computationally heavy problem per time step (having the most degrees of freedom) -even with an adaptive mesh.The higher densities have simpler and more localized electron dynamics, which demands a fine mesh with high resolution in only a small region of phase space.
In each case, simulations were performed for ten laser periods once the front of the laser pulse had reached the vacuum-plasma interface.For the adaptive solver, we used the adaption parameters w x,1 = w px,1 = 0.5 • r i/2 /f max , w x,2 = w px,2 = 0.5•r 2i /f max , w f = 6.25•10 −6 •r i /f max and δ thres = 10 −8 .These result in a relatively strongly refined mesh, even in regions with small values of the distribution function.Furthermore, the mesh was adapted every 20 time steps.In all cases, the code was compiled using the Intel compiler with -O3 optimizations and ran on a desktop PC using an Intel Xeon E3-1231 (3.4 GHz, 4 cores) CPU.
The performance of the adaptive solver can be quantified by the reduction of its memory footprint as well as the speedup compared to a uniform case.Figure 5(a) compares the memory footprints against system time, and Figure 5(b) shows the time it takes to perform a single time step, as a function of the time during the simulation (measured in laser periods).The reduction in memory footprint is time varying, but is at least a factor of 3 in all of the adaptive cases compared to the non-adaptive case.In all cases we observe a speedup of at least 4x (worst case, however for 50% of the simulation time we see ∼ 7x). the time it takes to advance the solution a time step as a function of time in laser periods (T ), when using an adaptive mesh to simulate a circularly polarized laser pulse impinging on plasmas with densities ne = 0.3nc (underdense plasma), ne = 1.5nc (SIT) and ne = 2.6nc (Hole Boring), as well as when using a uniform (Non-Adaptive) mesh.
The difference in time per time step is caused by a larger number of degrees of freedom in the lower density cases, however this difference is modest.
To compliment the results in Figure 5, we make a more in depth comparison of run time for the hole-boring case n e = 2.6n c , with and without adaption.The total runtime and division of the computational work between different parts of the solver is reported in Table 1.The runtime has been divided into four categories: (1) a time step category which includes the time for advancing the distribution function, (2) a category for the calculation of charge and current densities, (3) a category for the regridding procedure and (4) a category for the calculation of the electrostatic potential.The greatest proportion of the runtime, in all cases, was spent in one of these categories.
In both the uniform and adaptive cases, most of the run time is spent either on time stepping or calculating the charge and current.The contribution to the run time Table 1.A summary of the computational costs in key parts of veritas, with and without adaptive mesh refinement.Here, we limit the meaning of time step to only include the contribution from time advancement of the distribution functions and exclude any electromagnetic field or potential updates.

Category
Uniform from the solver for the transverse electromagnetic fields is marginal compared the time spent in the above regions.Furthermore, we notice that the work load distribution is mostly unaffected by the use of adaption.Notice that the computational cost associated with regridding was modest (< 2 %).However, there is a tradeoff between choosing small adaption parameters, resulting in a large number of points in the mesh, and high frequency for adaption, resulting in large computational costs for regridding.For example, in Ref. [22] a higher rate for adaption is used, with the consequence that a larger amount of the runtime is spent in regridding procedures.Optimal choices for the adaption parameters and regridding rate is an area to investigate to further optimize the accuracy and speed of the adaptive code.

Conclusions
The study of the Vlasov-Maxwell system of equations has long been of interest for modelling collisionless plasma dynamics.Eulerian methods are suitable for cases where high resolution of the distribution function is needed, but have so far been restricted by the computational expense.Here, an open source adaptive solver for the relativistic Vlasov-Maxwell system in one dimension, veritas, is presented.We have successfully demonstrated the solvers capabilities, leading to a speed-up of a factor of 7x in a typical scenario.It is shown that the adaptive approach is well suited for problems where a small component of the electron or ion populations are accelerated to high velocity.
The discretization of the Vlasov-equation in veritas is based on a high-order finite volume method, implementing the flux corrected transport algorithm to limit spurious oscillations in the distribution functions, in the presence of steep gradients.The use of efficient limiters is of critical importance to obtain a stable numerical scheme over a wide parameter range for self-consistent simulations with ultra intense fields, such as laser-plasma interactions.
The reduction of runtime and memory footprint using adaption comes from the reduction in the total number of cells in which the computation is performed, although with some additional overheads related to more complicated mesh structure and calculation of charge and current densities.In the cases and parameter sets considered in this paper, the overheads due to adaption were found to be minor.However, performance of the adaptive solver is also connected to the optimization of adaption parameters, e.g., the efficiency of the covering of marked cells using rectangles, and the threshold for as well as form of the refinement indicator.The optimal choices of these parameters are problem dependent.
Note that the local character of the explicit discretization in veritas makes it well suited for parallelization, resulting in an efficient tool for simulations of moderate laser intensities (with a 0 ∼ 1).However, in simulations of interaction of a plasma with ultra-relativistic fields (with a 0 1), performance for explicit methods is strongly limited by CFL-restrictions on the time step.Extensions to higher field strengths may therefore benefit from a more robust approach such as an explicit local time-stepping method [56] or the use of semi-Lagrangian methods which put fewer restrictions on the timestep [57].A third option may be an implicit treatment of the momentum derivatives in the Vlasov-equation, which would remove the severe scaling of the CFL-restriction due to the field strength; although this method will pose challenges for a parallel implementation.The small number of degrees of freedom in the adaptive mesh could be a significant advantage for the implicit approach.
Finally, concerning the extention of veritas from a 1D1P code to 2D2P: all numerical methods presented here are equally applicable.Based on the performance of the 1D1P code, we are optimistic about the potential performance gains that can be achieved by the use of an adaptive mesh in 2D2P.On the other hand, some further investigation is still needed with regard to the scaling of overheads and effect of managing the more complicated mesh structure in a distributed memory context.

Figure 2 .
Figure 2. Laser pulse interacting with an electron plasma slab, which is represented using a five level adaptive mesh.The color map shows the 10 largest orders of magnitude for the distribution function on a log 10 scale.The rectangular patches in each level are indicated by colored boxes.The coarsest level is not highlighted via a rectangle, but extends over the entire figure.The levels are nested, with each level (with yellow being the finest) contained within the next coarser level.

Figure 3 .
Figure3.Comparison of analytical solutions for density (dashed lines) and squared normalized vector potential (solid lines) and solutions that are calculated using veritas (black markers), for a circularly polarized laser pulse that impinges on an overdense plasma.The density has been normalized by 2nc, to fit on the same scale as the squared vector potential.The analytical results for different intensities are labelled with the colours red (a0 = 0.25), blue (a0 = 0.50), green (a0 = 0.75) and purple (a0 = 1.00).

Figure 4 .
Figure 4.The distribution function for electrons, calculated by veritas, at time t = 15T after the front of the laser pulse has reached the plasma, for the densities (a) ne = 0.3nc (in the underdense plasma regime), (b) ne = 1.5nc (in the SIT regime) and (c) ne = 2.6nc (in the hole-boring regime).The corresponding distribution function (grey) from PIC simulations is also shown.

n e = 1 .
5n c (in the SIT regime) and n e = 2.6n c (in the hole-boring regime).The n e = 0.3n c case is in the regime where the relativistic Raman and modulational instabilities merge[53,54] as evidenced by particle trapping and acceleration, the n e = 1.5n c case develops electron vortices on the front side which is in agreement with dynamics of the SIT regime[43,47,55] and the n e = 2.6n c is identified with characteristics of the hole-boring regime.The Figure also shows the distribution function calculated by the PIC code (grey) at the same instant of time.We observe remarkable agreement between epoch and veritas simulations.In the veritas simulations presented above, the dimensions of the simulation box were [0, 20]λ × [−8, 8]m e c for electrons and [0, 20]λ × [−400, 400]m e c for ions.Furthermore, the mesh consisted of five levels with refinement ratio r = 2 and n x = 292, n p = 150 points in the spatial and momentum directions respectively, on the coarsest level.The temperatures of both electrons and ions were taken to be T e = T i = 5 • 10 −4 m e c 2 .For epoch, a grid resolution of 200 cells per wavelength was used, with each cell spawning 8000 third order B-spline particles.

Figure 5 .
Figure 5. (a) Memory footprint as a function of system time and (b) the time it takes to advance the solution a time step as a function of time in laser periods (T ), when using an adaptive mesh to simulate a circularly polarized laser pulse impinging on plasmas with densities ne = 0.3nc (underdense plasma), ne = 1.5nc (SIT) and ne = 2.6nc (Hole Boring), as well as when using a uniform (Non-Adaptive) mesh.
At a coarse-fine interface, we define fluxes F fine,1 , F fine,2 and Fcoarse, where Fcoarse is inferred from F fine,1 , F fine,2 in order to ensure particle conservation.