Martin Campos Pinto

Towards smooth particle methods without smoothing

Keywords: Particle methods, Transport equations, A priori error estimates, Semi-Lagrangian methods, Remapped particle methods, Transport error indicators, Dynamic remapping strategies Mathematics Subject Classification (2000) 76M28, 35F10, 65M12

published or not. The documents may come

Towards smooth particle methods without smoothing

Introduction

Efficient and simple particle methods are a very popular tool for the numerical simulation of transport equations involved in many physical problems, ranging from fluid dynamics [START_REF] Chorin | Numerical study of slightly viscous flow[END_REF][START_REF] Cottet | Vortex Methods: Theory and Practice[END_REF] to kinetic (e.g., Vlasov) equations [START_REF] Hockney | Computer simulation using particles[END_REF][START_REF] Langdon | Plasma Physics via Computer Simulation[END_REF]. However, particle methods also suffer from weak convergence properties that cause difficulties in many practical cases. Specifically, it is known that they only converge in a strong sense when the particles present an extended overlapping, that is, when the number of overlapping particles tends to infinity as the mesh size h of their initialization grid tends to 0, see e.g. [START_REF] Beale | Vortex methods. II. Higher order accuracy in two and three dimensions[END_REF][START_REF] Raviart | An analysis of particle methods[END_REF]. Moreover, convergence rates are known to be suboptimal and to require vanishing moments for the particle shape functions, which prevents high orders to be achieved with positive shapes. In practice, extended particle overlapping is expensive and it involves an additional parameter to be optimized, such as the exponent q < 1 for which the particles radius behaves like h q ≫ h. In Particle-In-Cell (PIC) codes for instance, taking q < 1 typically leads to increasing the number of particles per cell faster than the number of cells, since the latter determines the effective radius of the particles [START_REF] Langdon | Plasma Physics via Computer Simulation[END_REF]. In Smoothed Particle Hydrodynamics (SPH) schemes it amounts to increasing the number of neighbors, i.e., interacting particles [START_REF] Rasio | Particle methods in astrophysical fluid dynamics[END_REF]. Because of these issues, accurate results often require numerically intensive runs, and in many cases the simulations do not meet the conditions of strong convergence. Significant oscillations are then produced, which are sometimes seen as a statistical noise that hampers interpretation of results, and can further cause large scale errors.

To suppress noise without resorting to extended particle overlapping, many methods (like the redeposition scheme introduced by Denavit [START_REF] Denavit | Numerical Simulation of Plasmas with Periodic Smoothing in Phase Space[END_REF] for plasma simulation and recently revisited as a Forward Semi-Lagrangian scheme (FSL), see e.g. [START_REF] Nair | A forward-trajectory global semi-Lagrangian transport scheme[END_REF][START_REF] Cotter | The remapped particle-mesh semi-Lagrangian advection scheme[END_REF][START_REF] Crouseilles | A forward semi-Lagrangian method for the numerical solution of the Vlasov equation[END_REF]) use periodic remappings, i.e., particle re-initializations that smooth out the evolution. However, frequent remappings sometimes introduce numerical dissipation which contradicts the benefit of using non-dissipative particle schemes, and to reduce this effect many works have been devoted to the design of accurate, low dissipative remapping schemes, including well-designed kernels and adaptive multi-resolution techniques, see e.g. [START_REF] Koumoutsakos | Inviscid Axisymmetrization of an Elliptical Vortex[END_REF][START_REF] Bergdorf | Multilevel Adaptive Particle Methods for Convection-Diffusion Equations[END_REF][START_REF] Bergdorf | A Lagrangian Particle-Wavelet Method[END_REF][START_REF] Wang | A Particle-In-Cell method with adaptive phase-space remapping for kinetic plasmas[END_REF][START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical flows at high Reynolds numbers[END_REF][START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF]. Alternative approaches have also been studied to improve the accuracy of the density evaluation without remapping the particles. In Beale's method [START_REF] Beale | On the accuracy of vortex methods at large times[END_REF] for instance, new weights are iteratively computed from the positions of the particles to evaluate the density by an approximate interpolation technique. In Cottet's scheme [START_REF] Cottet | Artificial Viscosity Models for Vortex and Particle Methods[END_REF] the particle weights are also modified over time and in order to reduce the resulting smoothing, this operation is performed by introducing the proper amount of anisotropic viscosity so as to minimize the dominant (antidiffusive) part of the truncation error. To avoid smoothing the solution in a non-reversible way, Strain's method [START_REF] Strain | 2D Vortex Methods and Singular Quadrature Rules[END_REF] does not modifies the particle weights but dynamically constructs local interpolation formulas using adaptive cells that contain a given number of neighboring particles. In [START_REF] Cohen | Optimal Approximations of Transport Equations by Particle and Pseudoparticle Methods[END_REF] this idea is further investigated to establish optimal error estimates. Some authors have also studied the effect of using less regular convolution kernels, see e.g. [START_REF] Holm | Euler-alpha and vortex blob regularization of vortex filament and vortex sheet motion[END_REF].

In this article we present a new class of particle schemes where polynomial expansions of the characteristic flow are locally computed to transport the numerical particles with improved accuracy and less remappings. Specifically, our method modifies the particle shape functions using polynomial transformations which coefficients involve local derivatives of the backward flow: the first order version is a linearly-transformed particle (LTP) method, the second order case is a quadratically-transformed particle (QTP) scheme, and so on. With this approach the particles do not present an extended overlapping, and apart from the occasional remapping steps their weights are not modified over time. Moreover, the method is flexible in the sense that several "kernels" can be used for the particle shape functions.

Our contributions are twofold. On a practical level we provide explicit implementations of the LTP and QTP schemes that comply with the algorithmic structure of most existing particle codes, in the sense that they only involve pointwise evaluations of the forward characteristic flow. We show that the resulting schemes are not subject to a CFL-like condition and also propose a dynamic criterion to automatically select the time steps where the particles should be remapped.

On a theoretical level we extend these particle methods up to arbitrary polynomial orders and prove that they converge in the uniform norm at a rate that depends on the smoothness of the characteristic flow, provided that the initial particle shape functions are Lipschitz. Our proof does not require remappings, extended overlapping or vanishing moments for the particles.

Finally, numerical tests using different transport problems demonstrate the accuracy of the proposed methods compared to standard particle schemes, and their robustness with respect to the remapping periods. In particular, it is shown that QTP particles can be transported without remappings on very long periods of times, without hampering the accuracy of the numerical solutions.

Deforming the particles is not a new idea. Our method can be viewed as a modified version of Hou's method [START_REF] Hou | Convergence of a Variable Blob Vortex Method for the Euler and Navier-Stokes Equations[END_REF] further developed in [START_REF] Cottet | Vortex Methods with Spatially Varying Cores[END_REF][START_REF] Bergdorf | Multilevel Adaptive Particle Methods for Convection-Diffusion Equations[END_REF], where particles are deformed using a global deformation mapping rather than local approximations of the backward flow. Variants have been studied by other authors. In [START_REF] Cohen | Optimal Approximations of Transport Equations by Particle and Pseudoparticle Methods[END_REF] Cohen and Perthame observed that by linearizing the flow around the particle trajectories one obtains a convergent method (in L 1 ) with particles scaled with their initialization grid, and no remappings. More recently, Alard and Colombi [START_REF] Alard | A cloudy Vlasov solution[END_REF] derived a "cloudy" method similar to ours by evolving Gaussian particles with locally affine force fields in PIC simulations of the Vlasov-Poisson system.

The outline is as follows. After a rapid overview of the main particle methods in Section 2 where we introduce some notations, we describe explicit schemes in Section 3: In Section 3.1 we recall classical approximation algorithms used to initialize and remap the particles on a cartesian grid, and in Sections 3.2 and 3.4 we provide discrete schemes for the LTP and QTP schemes that only involve pointwise evaluations of the forward flow. In Section 3.5 we present an original criterion to dynamically select the remapping time steps, based on the theoretical study developed in Section 4. There we extend our particle methods to arbitrary orders, and establish a priori error estimates. Numerical results are eventually presented in Section 5, that validate the robustness of our approach.

A brief review of particle methods

Following [START_REF] Raviart | An analysis of particle methods[END_REF][START_REF] Cohen | Optimal Approximations of Transport Equations by Particle and Pseudoparticle Methods[END_REF] we consider the transport equation

∂ t f (t, x) + u(t, x) • ∇f (t, x) = 0, t ∈ [0, T ], x ∈ R d (1)
associated with an initial data f 0 : R d → R, a final time T and a velocity field u : [0, T ]×R d → R d . In fluid problems for instance we have d = 2, 3, while in kinetic formulations R d is a phase space with d ≤ 6 and u is a generalized velocity field with components of velocity and acceleration. We assume that u is smooth enough (e.g., Lipschitz) for the characteristic trajectories X(t) = X(t; s, x), solutions to

X ′ (t) = u(t, X(t)), X(s) = x, (2)
to be defined on [0, T ] for all x ∈ R d and s ∈ [0, T ]. In particular, the characteristic flow F s,t : x → X(t) is invertible and satisfies (F s,t ) -1 = F t,s . Solutions to (1) read then

f (t, x) = f 0 ((F 0,t ) -1 (x)) for t ∈ [0, T ], x ∈ R d . (3)
Let the time steps be denoted by t n = n∆t for some ∆t > 0. Throughout this article we make the simplifying assumption that the exact forward flow is known,

F n ex := F t n ,t n+1 : R d → R d (4)
or some accurate approximation to it that we denote by F n , see in particular Section 3. To apply our method to non-linear problems one must then provide a scheme that computes F n from the particle approximation to f (t n ). In [START_REF] Pinto | Noiseless Vlasov-Poisson simulations with linearly transformed particles[END_REF] this has been done for the 1d1v Vlasov-Poisson equation.

For simplicity, we also restrict ourselves to incompressible fields satisfying div u = 0. In this case the flow is measure preserving in the sense that its Jacobian matrix

J F s,t (x) = ∂ j (F s,t ) i 1≤i,j≤d satisfies det J F s,t (x) = 1 for s, t ∈ [0, T ], x ∈ R d . (5)
2.1 The traditional smoothed particle method (TSP)

In the standard "academic" particle method [START_REF] Raviart | An analysis of particle methods[END_REF], numerical solutions are typically computed as follows: considering deterministic initializations for simplicity, the initial data f 0 is first approximated by a collection of particles on a regular (say, cartesian) grid of step h > 0,

f 0 h,ǫ (x) := k∈Z d w k (f 0 )ϕ ǫ (x -x 0 k ) with x 0 k := hk
and with weights typically defined as

w k (f 0 ) := x 0 k +[-h 2 , h 2 ] d f 0 (x) dx or w k (f 0 ) := h d f 0 (x 0 k ). (6)
Here ϕ ǫ = ǫ -d ϕ(•/ǫ) is a particle shape function with radius proportional to ǫ, usually seen as a smooth approximation of the Dirac measure obtained by scaling a compactly supported "cut-off" function ϕ for which common choices include B-splines and smoothing kernels with vanishing moments, see e.g. [START_REF] Koumoutsakos | Inviscid Axisymmetrization of an Elliptical Vortex[END_REF][START_REF] Cottet | Vortex Methods: Theory and Practice[END_REF]. Particle centers are then pushed forward by following the numerical flow F n , leading to

f n+1 h,ǫ (x) := k∈Z d w k (f 0 )ϕ ǫ (x -x n+1 k ) ≈ f (t n+1 , x) with x n+1 k := F n (x n k ).
In the classical error analysis [START_REF] Beale | Vortex methods. II. Higher order accuracy in two and three dimensions[END_REF][START_REF] Raviart | An analysis of particle methods[END_REF], the above process is seen as (i) an approximation (in the distribution sense) of the f 0 by a collection of weighted Dirac measures, (ii) the exact transport of the Dirac particles, and (iii) the smoothing of the resulting k w k (f 0 )δ x n k with the convolution kernel ϕ ǫ . The classical estimate reads then as follows: if for some prescribed integers m > 0 and r > 0, the cut-off ϕ has m orders of smoothness and satisfies a moment condition of order r, namely if ϕ = 1, |y| r |ϕ(y)| dy < ∞ and

y s 1 1 . . . y s d d ϕ(y 1 , . . . , y d ) dy = 0 for s ∈ N d with 1 ≤ s 1 + • • • + s d ≤ r -1,
then there exists a constant C independent of f 0 , h or ǫ, such that

f (t n ) -f n h,ǫ L µ ≤ C ǫ r f 0 W r,µ + (h/ǫ) m f 0 W m,µ , 1 ≤ µ ≤ ∞. (7)
Recently, Cohen and Perthame [START_REF] Cohen | Optimal Approximations of Transport Equations by Particle and Pseudoparticle Methods[END_REF] observed that using weights defined by

w k (f 0 ) := R d f 0 (x) φh (x -x 0 k ) dx (8)
where φh = φ(•/h) is derived from a compactly supported φ satisfying

k∈Z d k s 1 1 . . . k s d d φ(y -k) = y s 1 1 . . . y s d d for s ∈ N d with 0 ≤ s 1 + • • • + s d ≤ m -1,
one has the improved estimate

f (t n ) -f n h,ǫ L µ ≤ C ǫ r f 0 W r,µ + (h/ǫ) m f 0 L µ (9)
with a constant that is again independent of f 0 , h or ǫ. Note that ( 9) improves (7) since m is not constrained by the smoothness of f 0 , which allows to reach higher convergence rates. Indeed balancing the error terms in the above estimates suggests to take ǫ ∼ h q with q = m m+r , yielding a convergence in h q = h rm m+r . In particular, if f 0 ∈ W s,µ for some integer s then the best possible rate with standard weights is only h s/2 f 0 W s,µ , obtained with m = r = s. With the improved weights instead, one can take a higher value for m and obtain estimates close to h s f 0 W s,µ . Moreover, the latter approach also allows to improve (i.e., reduce) the particle overlapping, since the corresponding exponents are q = 1 2 and m m+r ≈ 1, respectively.

In either case, we see from the terms ǫ r ∼ h qr in the estimates that extended particle overlapping does not only make the simulations more expensive, it also deteriorates their convergence order.

The forward semi-Lagrangian scheme (FSL)

In forward semi-Lagrangian schemes (also called remapped or remeshed particle methods), particles typically have the scale of their initialization grid (ǫ = h) but they are periodically remapped on that grid, say every N r time steps. Thus, if

A h : g → k∈Z d w k (g)ϕ h (x -x 0 k ), x 0 k := hk
is an approximation operator on the grid, and

T n h : ϕ h (• -x n k ) → ϕ h (• -F n (x n k
)) is the fixed-shape particle transport operator, FSL schemes take the generic form [START_REF] Chui | A characterization of multivariate quasi-interpolation formulas and its applications[END_REF] In practice these schemes achieve high levels of accuracy when the remapping operators are well designed, see e.g. [START_REF] Bergdorf | A Lagrangian Particle-Wavelet Method[END_REF][START_REF] Wang | A Particle-In-Cell method with adaptive phase-space remapping for kinetic plasmas[END_REF][START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical flows at high Reynolds numbers[END_REF][START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF] for recent applications to challenging problems in plasmas physics and fluid mechanics.

f n+1 h = k∈Z d wn k ϕ h (• -F n (x n k )) = T n f n h with f n h := A h f n h if n ∈ N r N f n h otherwise.

Particle transport with polynomial transformations

In this article we develop a lesser-known approach where particles are subject to polynomial transformations to better approximate the backward flow involved in the exact transport [START_REF] Beale | On the accuracy of vortex methods at large times[END_REF]. At the first order for instance, this leads to formally defining linearly-transformed particles as

ϕ n h,k (x) := ϕ h D n k (x -x n k ) . (11)
Thus, in addition to pushing forward the particle centers one needs to compute d × d deformation matrices D n k , k ∈ Z d , that approximate the local Jacobian matrices of the backward flow. Similarly, at the second order the particle shape functions are transformed with local quadratic mappings which coefficients involve the second derivatives of the backward flow, and so on.

Because in practice it is often necessary to remap the particles to obtain satisfactory results, we shall write our particle methods in a form similar to [START_REF] Chui | A characterization of multivariate quasi-interpolation formulas and its applications[END_REF]. However, it will be shown in Section 4 that remappings are not required for the convergence of the method. An important consequence is that the optimal remapping frequency is significantly reduced compared to fixed-shape particle schemes, and the method is free of any CFL-like condition.

In the sequel we let x ∞ := max i |x i | for vectors and A ∞ := max i j |A i,j | for matrices. For v in the Sobolev space W m,∞ (ω) with index m > 0 we write

|v| m,ω := max i d l 1 =1 • • • d l m =1 ∂ l 1 • • • ∂ l m v i L ∞ (ω) (12)
and for conciseness we shall drop the domain ω ⊂ R d when it is the whole space.

3 Particle transport with linear and quadratic shape transformations

In this section we describe numerical schemes to transport a collection of particles

f 0 h = k∈Z d w k ϕ 0 h,k (13)
with particle-wise operators approximating the exact transport (3), denoted here

T n ex : ϕ 0 h,k → ϕ 0 h,k ( Bn ex (•)) where Bn ex := ( F n ex ) -1 , F n ex := F 0,t n . ( 14
)
Here and below we use a bar to distinguish flows and transport operators defined on the global time interval [0, t n ] from those on the single time step [t n , t n+1 ], see (4). To comply with the algorithmic structure of standard particle codes, our schemes only involve pointwise evaluations of a given forward flow F n which may either be the exact flow [START_REF] Bergdorf | Multilevel Adaptive Particle Methods for Convection-Diffusion Equations[END_REF] or more likely some approximation to it. In practice F n will typically be given by the algorithm used to push the particles in an existing code.

For simplicity we consider particles initially centered on a cartesian grid with nodes x 0 k := hk, k ∈ Z d , scaled to the grid and normalized in L 1 , i.e.,

ϕ 0 h,k (x) := ϕ h (x -x 0 k ) := h -d ϕ(h -1 x -k). (15)
The crucial point is the introduction of a deformation matrix D n k which is a d × d matrix representing the local Jacobian of the (backward) flow. And in the case of quadratic transformations, additional matrices (Q n k ) i , 1 ≤ i ≤ d, are introduced to represent the local Hessian matrices of the flow. Here ϕ is a reference shape supported on the d-dimensional cube B ℓ ∞ (0, ρ 0 ) (see Section 3.1 below), so that the particles are initially supported in small cubes

Σ 0 h,k := supp(ϕ 0 h,k ) = x 0 k + supp(ϕ h ) ⊂ B ℓ ∞ (x 0 k , hρ 0 ) ( 16
)
and with the shape transformations these supports undergo polynomial deformations. For instance, with linearly-transformed particles their expression becomes

Σ n h,k := supp(ϕ n h,k ) = x n k + (D n k ) -1 (supp(ϕ h )),
see Section 3.2 below.

Particle approximations using a cartesian grid

Several choices can be made for the reference shape function ϕ in [START_REF] Cottet | Vortex Methods with Spatially Varying Cores[END_REF]. One standard option consists of taking the interpolating kernel M ′ 4 of Monaghan [START_REF] Monaghan | Extrapolating B. Splines for Interpolation[END_REF],

ϕ(x) = d i=1 M ′ 4 (x i ) with M ′ 4 (x i ) =      1 -5 |x i | 2 2 + 3 |x i | 3 2 if 0 ≤ |x i | ≤ 1 1 2 (2 -|x i |) 2 (1 -|x i |) if 1 ≤ |x i | ≤ 2 0 otherwise. ( 17

) Particles are then initialized and remapped by a call to the approximation operator

A h : g → k∈Z d w k (g)ϕ 0 h,k with w k (g) := h d g(x 0 k ). ( 18
)
The resulting A h is an interpolation that reproduces second-degree polynomials, hence it is third order accurate.

Another option is to use cardinal B-splines, recursively defined for x ∈ R by

B 0 (x) := 1 -1 2 , 1 2 (x) and B p (x) := (B p-1 * B 0 )(x) = x+ 1 2 x- 1 2 B p-1 , p ≥ 1, so that B 1 (x) = max{1-|x|, 0}
is the traditional "hat-function", B 3 is the centered cubic B-spline, and so on. For the reference particle shape function we then set

ϕ(x) := d i=1 B p (x i ) supported on supp(ϕ) = B ℓ ∞ (0, ρ 0 ) with ρ 0 := p+1 2 .
(19) For the initialization and remappings we can then use standard approximation schemes that rely on the fact that the span of their integer translates (15) contains the space P p of polynomials with coordinate degree less or equal to p. Specifically, we can use the quasi-interpolation schemes described by [START_REF] Chui | A characterization of multivariate quasi-interpolation formulas and its applications[END_REF] and [START_REF] Unser | On the approximation power of convolution-based least squares versus interpolation[END_REF], where highorder B-spline approximations are locally obtained by pointwise evaluations of the target function. The resulting approximation reads

A h : g → k∈Z d w k (g)ϕ 0 h,k with w k (g) := h d l ∞ ≤m p a l g(x 0 k+l ), a l := 1≤i≤d a l i
with symmetric coefficients a l = a -l computed with the iterative algorithm in [10, Section 6]: for the first odd orders we obtain

• m p = 0 and a 0 = 1 for p = 1,

• m p = 1 and (a 0 , a 1 ) = ( 86 , -1 6 ) for p = 3, • m p = 4 and (a 0 , a 1 , a 2 , a 3 , a 4 ) = ( 503288 , -1469 3600 , 7 225 , 13 3600 , 1 14400 ) for p = 5. The resulting A h reproduces polynomials in P p , hence it is of order p + 1: we have

A h g -g L ∞ ≤ c A h q+1 |g| q+1 for q ≤ p (20)
with a constant c A independent of h.

In Section 5 we will show numerical results obtained by applying our transport schemes to either M ′ 4 particles ( 17), or cubic spline particles [START_REF] Holm | Euler-alpha and vortex blob regularization of vortex filament and vortex sheet motion[END_REF] with p = 3. We already note that the latter choice results in remappings that are fourth-order, but are also more numerically dissipative due to the wider stencil.

Particle transport with linear shape transformations

The linearly-transformed particle (LTP) method essentially consists of transporting the particles along local linearizations of the exact characteristic flow. Specifically, it is based on the observation that the operator

T n h,(1) : ϕ h (• -x 0 k ) → ϕ h D n k (• -x n k ) with x n k = F n ex (x 0 k ) D n k = J Bn ex (x n k ) (21)
corresponds to the exact transport of the k-th particle along the affine flow

F n (1),k (x) := x n k + (D n k ) -1 (x -x 0 k ). ( 22
)
Indeed the resulting particles satisfy

ϕ h D n k (x -x n k ) = ϕ h ( Bn (1),k (x) -x 0 k ) with Bn (1),k (x) := ( F n (1),k ) -1 (x) = x 0 k + D n k (x -x n k ), (23)
and by deriving 5), we verify that ( 22) is the linearization of the forward flow F n ex at x 0 k . The approximate density reads then

J Bn ex ( F n ex (x 0 k )) = (J F n ex (x 0 k )) -1 from (
f n h (x) = k∈Z d w k ϕ n h,k (x) with ϕ n h,k (x) = ϕ h D n k (x -x n k ) . (24)
In practice we propose to implement approximate versions of ( 21), i.e.,

T n h,ltp : ϕ 0 h (• -x 0 k ) → ϕ h D n k (• -x n k ) with x n k ≈ F n ex (x 0 k ) D n k ≈ J Bn ex (x n k ), (25)
that only involve pointwise evaluations of F n ≈ F n ex . For the centers we simply set

x n k := F n-1 (x n-1 k ) = F n (x 0 k ) where F n := F n-1 • • • • • F 0 .
For the deformation matrices D n k we can think of different schemes, the benefits of which may depend on the context. The first one may be called direct as it consists of directly computing them when needed, using the position of the neighboring particles. Specifically, we approximate the derivatives of F n ≈ F n ex using finite differences on the original cartesian grid. With a centered formula this gives

Jn k := (x n k+e j -x n k-e j ) i 2h 1≤i,j≤d ≈ J F n ex (x 0 k ) (26)
and using the formal identity

J F n (x 0 k )J Bn (x n k ) = I d , we set D n k := ( Jn k ) -1 . (27)
The second scheme is incremental: the deformation matrices are carried with the particles and updated at each time step with a few calls to F n . Using the formal identities

J Bn+1 (x n+1 k ) = J Bn (x n k )J B n (x n+1 k ) = J Bn (x n k )J F n (x n k ) we compute J n k := F n i (x n k + he j ) -F n i (x n k -he j ) 2h 1≤i,j≤d ≈ J F n ex (x n k )
where e j := (δ j,l ) 1≤l≤d , and then set

D n+1 k := D n k (J n k ) -1 (28)
see Remark 4 below for an example in the case of a simple Vlasov-Poisson equation. Note that since B0 ex = I, the deformation matrices are initialized with D 0 k := I d .

Remark 1 A rigorous study of the conditions under which the matrices Jn k or J n k are invertible is not included here, but we note that since they approximate Jacobian matrices with determinant equal to 1 they are not likely to be singular, and in the numerical tests presented in Section 5 this situation never occurred.

A numerical indicator for the local flow error

Before turning to second order shape transformations, we observe that it is possible to define numerical indicators for the local flow error

e n B,k := Bn (1),k -Bn ex L ∞ (Σ n h,k )
where

Σ n h,k = supp(ϕ n h,k ). ( 29
)
Indeed, if we assume that the numerical flow F n is an accurate approximation to F n ex , then we have Bn ex (x n k+ℓ ) = Bn ex ( F n (x 0 k+ℓ )) ≈ x 0 k+ℓ with good accuracy, so that we may estimate the local flow error (29) by ên B,(1),k := max

ℓ ∞ ≤1 Bn (1),k (x n k+ℓ ) -x 0 k+ℓ ∞ . ( 30
)
Note that when the matrix D n k is computed with the scheme ( 26)-( 27), it is possible to express the error Bn

(1),k (x n k+ℓ )-x 0 k+ℓ in terms of finite difference approximations to the second order derivatives of F n . Indeed, the definition of Bn

(1),k yields

Bn (1),k (x n k+ℓ ) -x 0 k+ℓ = D n k (x n k+ℓ -x n k ) -(x 0 k+ℓ -x 0 k ) = D n k x n k+ℓ -x n k -Jn k hℓ = D n k x n k+ℓ -x n k - d j=1 ℓ j 2 (x n k+e j -x n k-e j ) .
In particular, for ℓ = ±e i (both values give the same result) we obtain

Bn (1),k (x n k+ℓ ) -x 0 k+ℓ = h 2 2 D n k x n k+e i -2x n k + x n k-e i h 2
which involves a finite difference approximation of the second derivative

∂ 2 i F n . For ℓ = e i + e l , we find Bn (1),k (x n k+ℓ ) -x 0 k+ℓ = h 2 2 D n k (x n k+e i +e l -x n k+e i ) -(x n k -x n k-e l ) h 2 + (x n k+e i +e l -x n k+e l ) -(x n k -x n k-e i ) h 2
which involves an approximation of the second derivative (∂ i + ∂ l ) 2 F n , and so on. Finally, from the quantities [START_REF] Unser | On the approximation power of convolution-based least squares versus interpolation[END_REF] we derive a global indicator for the flow error, ên B,(1) := sup

k∈Z d ên B,(1),k . (31)

Particle transport with quadratic shape transformations

Extending the above idea we define a quadratically-transformed particle (QTP) method based on second-order expansions of the exact flow around the particle trajectories. Because the exact transport operator ( 14) involves the backward flow Bn ex , we choose to define a quadratic approximation of the latter at x n k , denoted

Bn (2),k (x) := x 0 k + D n k (x -x n k ) + 1 2 (x -x n k ) t Q n k (x -x n k ). ( 32
)
Here the particle centers and deformation matrices may be computed as above.

As for the quadratic deformation terms

(x -x n k ) t Q n k (x -x n k ) = d j 1 ,j 2 =1 (Q n k ) i (x -x n k ) j 1 (x -x n k ) j 2 1≤i≤d
, they involve approximations of the Hessian matrices of the backward flow,

(Q n k ) i ≈ H ( Bn ex ) i (x n k ) := ∂ j 1 ,j 2 ( Bn ex ) i (x n k ) 1≤j 1 ,j 2 ≤d ( 33
)
and to compute them in practice we shall give again two discrete schemes that only involve pointwise evaluations of the numerical forward flow F n .

To define the resulting particles however, we observe that care must be taken. Indeed, the flow (32) being quadratic it is not invertible and hence there is no guarantee that the support of ϕ h ( Bn

(2),k (•) -x 0 k ) is contained in a ball of radius O(h).
In practice this causes strong oscillations in the numerical density, which follows from the fact that the expansion ( 32) is only relevant in a small region close to x n k . For this reason it is necessary to restrict the numerical particles on carefully shaped domains Σ n h,k ≈ F n ex (Σ 0 h,k ), see [START_REF] Crouseilles | A forward semi-Lagrangian method for the numerical solution of the Vlasov equation[END_REF]. In Theorem 2 below we will show that for the second-order convergence of the method one may take

Σ n h,k := F n (1),k (B ℓ ∞ (x 0 k , hρ n h,k )) with ρn h,k := ρ 0 + 1 h Bn (1),k -Bn ex L ∞ ( F n ex (Σ 0 h,k ))
and to estimate the terms ρn h,k we can use the local indicators [START_REF] Unser | On the approximation power of convolution-based least squares versus interpolation[END_REF]. In practice however, we have observed that these supports were sometimes too large, yielding visible oscillations in the solutions. And because the latter are likely caused by the non-invertibility of the quadratic flow Bn

(2),k , we decided to further restrict the particle supports to the regions where this flow is locally invertible. Specifically, we define the a-priori particle supports as

Σ n h,k := {x ∈ F n (1),k (B ℓ ∞ (x 0 k , hρ n h,k )) : det(J Bn (2),k (x)) > 0}. Since J Bn (2),k (x) = D n k + j ′ (Q n k ) i,j,j ′ (x -x n k ) j ′ 1≤i
,j≤d this strategy is easy to implement and in Section 5 we will see that it results in a very robust numerical method. In summary, in the QTP method the numerical density is defined as

f n h (x) = k∈Z d w k ϕ n h,k (x) with ϕ n h,k (x) = 1 Σ n h,k (x)ϕ h D n k (x -x n k ) + 1 2 (x -x n k ) t Q n k (x -x n k ) , (34)
and it only remains to specify how the quadratic deformation matrices (33) are computed in practice. As in Section 3.2 we describe two implementations. In the direct approach we first compute approximate forward Hessian matrices with finite differences on the original grid, ( Hn

k ) i := (h) -2 1 α 1 ,α 2 =0 (-1) α 1 +α 2 x n k+α 1 e j 1 +α 2 e j 2 i 1≤j 1 ,j 2 ≤d ≈ H ( F n ex ) i (x 0 k ).
Differentiating twice (at x 0 k ) the formal identity I = Bn F n we then obtain

0 = J F n (x 0 k ) t H ( Bn ) i (x n k )J F n (x 0 k ) + d j=1 J Bn (x n k ) i,j H ( F n ) j (x 0 k ) (35)
so that we finally set

(Q n k ) i := -(D n k ) t d j=1 (D n k ) i,j ( Hn k ) j D n k . (36)
In the incremental approach the d × d matrices (Q n k ) i are stored with the particles. Since B0 ex = I they are initialized with (Q 0 k ) i := 0. To update them we differentiate twice (at x n+1 k ) the formal identity Bn+1 = Bn B n : this gives

H ( Bn+1 ) i (x n+1 k ) = J B n (x n+1 k ) t H ( Bn ) i (x n k )J B n (x n+1 k ) + d j=1 J Bn (x n k ) i,j H (B n ) j (x n+1 k ) (37)
and we observe that rewriting (35) on the local time step yields a formula expressing the local backward Hessian matrices in terms of the forward ones. Thus we first compute finite difference approximations of the matrices

H (F n ex ) i (x n k ), (H n k ) i := (h) -2 1 α 1 ,α 2 =0 (-1) α 1 +α 2 (F n ) i (x n k + h(α 1 e j 1 + α 2 e j 2 )) 1≤j 1 ,j 2 ≤d
then we approximate the local backward Hessian matrices similarly as in (36), ( Ȟn

k ) i := -( Jn k ) t d j=1 ( Jn k ) i,j (H n k ) j Jn k ≈ H (B n ex ) i (x n+1 k )
where Jn k = (J n k ) -1 is the approximate local backward Jacobian. Finally using (37) we let

(Q n+1 k ) i := ( Jn k ) t (Q n k ) i Jn k + d j=1 (D n k ) i,j ( Ȟn k ) j = ( Jn k ) t (Q n k ) i - d j,j ′ =1 (D n k ) i,j ( Jn k ) j,j ′ (H n k ) j ′ Jn k . (38)
3.5 A dynamic criterion to select the remapping steps An interesting by-product of the error analysis developed in Section 4 is that it is possible to propose a numerical criterion for automatically selecting the time steps where the particles should be remapped. Although our criterion does not always perform as well as the optimal static strategy, it gives very satisfactory results when tested on our different benchmark problems, and it is local in the sense that it only involves communications between neighboring particles.

To describe it let us denote by n i , i = 0, . . . , R -1, the time steps of the initial approximation and further remappings, as selected by a given strategy. Then if

T n,m h = T m-1 h • • • T n
h denotes the transport operator acting without remappings between the times t n and t m , the remapped particle scheme S N h that maps f 0 to f N h reads

S N h f 0 = f N h = f n R h = T n R-1 ,n R h A h f n R-1 h = (T n R-1 ,n R h A h ) • • • (T n 0 ,n 1 h A h )f 0
(where n R = N ) and the global error

e N h := S N h -T 0,N ex L ∞ satisfies e N h ≤ (T n R-1 ,n R h -T n R-1 ,n R ex )A h f n R-1 h L ∞ + T n R-1 ,n R ex (A h -I)f n R-1 h L ∞ + T n R-1 ,n R ex (S n R-1 h -T 0,n R-1 ex )f 0 L ∞ ≤ (T n R-1 ,n R h -T n R-1 ,n R ex )A h f n R-1 h L ∞ + (A h -I)f n R-1 h L ∞ + e n R-1 h ≤ R-1 i=0 (T n i ,n i+1 h -T n i ,n i+1 ex )A h f n i h L ∞ + (A h -I)f n i h L ∞ .
(here f 0 h = f 0 ). Thus e N h essentially consists of the transport and remapping errors. Our heuristic is then as follows: although [START_REF] Hou | Convergence of a Variable Blob Vortex Method for the Euler and Navier-Stokes Equations[END_REF] tells us that the remapping errors can grow quickly when the smoothness of the solutions deteriorate, in practice we have observed that they do not depend much on the selected remapping steps. As for the transport errors, their increase is comparatively fast when the remapping period grows large, and it is easily seen that by resetting the associated flow to the identity, each remapping resets the transport error to 0. Therefore it seems reasonable to remap the particles when the estimated transport error becomes larger than the estimated remapping error. Specifically we shall remap f n h when

C remap E (T n i ,n h -T n i ,n ex )A h f n i h ) ≥ E (A h -I)f n h ) . (39)
Here n i denotes the last remapping step preceding n and C remap is a parameter to be determined from numerical experiments. It remains to specify numerical indicators for the transport and remapping errors involved in the remapping criterion (39). To do so we rely on the estimates that will be derived in Section 4. From (61) and ( 69)-(71) we derive the indicator

E (T n i ,n h -T n i ,n ex )A h f n i h ) := 1 + ên B,(1) (h) h d ên B,(r) (h) h f n i h L ∞ (40)
with r = 1 in the LTP case and 2 in the QTP case. Here the indicators for the backward flow errors are computed as in [START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical flows at high Reynolds numbers[END_REF], using the flows ( 23) and [START_REF] Wang | A Particle-In-Cell method with adaptive phase-space remapping for kinetic plasmas[END_REF], respectively. As for the remapping error, we use a first order estimate for simplicity

(A h -I)f n h L ∞ h|f n h | 1 h d j=1 sup k∈Z d |∂ j f n h (x n k )| ( 41
)
and to approximate the spatial derivatives we write

∂ j f n h (x n k ) ≈ ∂ j (f n i h • B n i ,n ex )(x n k ) = d ℓ=1 ∂ ℓ f n i h (B n i ,n ex (x n k ))∂ j (B n i ,n ex ) ℓ (x n k ).
Observing next that for the remapped particle scheme the approximations involved in [START_REF] Monaghan | Extrapolating B. Splines for Interpolation[END_REF]

read x n k ≈ F n i ,n ex (x 0 k ) and D n k ≈ J B n i ,n ex (x n k )
, we derive the following indicator for the remapping error,

E (A h -I)f n h ) := h d j=1 sup k∈Z d d ℓ=1 ∂ ℓ f n i h (x 0 k )(D n k ) ℓ,j . (42)
To estimate the spatial derivatives of f n i h we finally suggest to use finite differences on the cartesian grid {x 0 k = hk : k ∈ Z d }, since the point values f n i h (x 0 k ) are computed in the remapping algorithm. The numerical performances of the resulting dynamic remapping strategy will be assessed in Section 5.3.

Error analysis and higher order particle transformations

In order to cast the linear [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF] and quadratic (34) particle methods into a general framework, we now consider a transport operator of the form

T n h ϕ 0 h,k (x) := 1 Σ n h,k (x)ϕ 0 h,k ( Bn h,k (x)). ( 43
)
Here Bn h,k ≈ Bn ex is the approximated backward flow for the k-th particle, 1 is the set characteristic function and Σ n h,k is an a priori support for the particle. As explained in Section 3.4, such supports are needed when the domains ( Bn h,k ) -1 (Σ 0 h,k ) are not easily computable, or when they are very large compared with F n ex (Σ 0 h,k ). For linearly-transformed particles this is not the case, and by a proper choice of Σ n h,k the transport operator (43) can be simplified into

T n h ϕ 0 h,k (x) := ϕ 0 h,k ( Bn h,k (x)). (44)

Two preliminary estimates

It will be useful to state preliminary estimates based on the overlapping constant

Θ n (h) := sup x∈R d # {k ∈ Z d : x ∈ Σ n h,k } (45)
and on the backward flow error

e n B (h) := sup k∈Z d Bn h,k -Bn ex L ∞ (Σ n h,k ) . (46)
We first prove an estimate for the general operator (43).

Lemma 1 If the exactly transported particles vanish outside the domains Σ n h,k ,

F n ex (Σ 0 h,k ) ⊂ Σ n h,k , k ∈ Z d , (47)
then the approximate transport operator (43) satisfies

( T n h -T n ex )f 0 h L ∞ e n B (h)Θ n (h) h f 0 L ∞
with a constant independent of h, f 0 , u and n.

Proof For x ∈ R d , we let K n h (x) := {k ∈ Z d : x ∈ Σ n h,k } and infer from (47) that

T n h ϕ 0 h,k (x) -ϕ 0 h,k ( Bn ex (x)) = 0 if k / ∈ K n h (x) ϕ 0 h,k ( Bn h,k (x)) -ϕ 0 h,k ( Bn ex (x)) otherwise.
It follows that

|( T n h -T n ex )f 0 h (x)| = k∈Z d w k T n h ϕ 0 h,k (x) -f 0 h ( Bn ex (x)) = k∈K n h (x) w k ϕ 0 h,k ( Bn h,k (x)) -ϕ 0 h,k ( Bn ex (x)) ≤ k∈K n h (x) |w k ||ϕ 0 h,k | 1 Bn h,k -Bn ex L ∞ (Σ n h,k ) ≤ C f 0 L ∞ h -1 Θ n (h)e n B (h),
where we have used the scaling |ϕ 0 h,k | 1 ∼ h -1-d and the fact that particle weights defined as in Section 3.1 satisfy

|w k | h d f 0 L ∞ for some absolute constant. ⊓ ⊔
If we next consider the simpler operator (44) and denote

Σ n h,k := F n ex (Σ 0 h,k ) ∪ F n h,k (Σ 0 h,k ), (48)
then the definitions (44) and (43) are equivalent. Moreover Assumption (47) is readily fulfilled and it is possible to establish an a priori bound for the corresponding overlapping constant (45) and the transport error, that is either based on the above backward flow error (46) or on the forward flow error

e n F (h) := sup k∈Z d F n h,k -F n ex L ∞ (Σ 0 h,k ) . (49)
Lemma 2 The approximate transport operator (44) satisfies

( T n h -T n ex )f 0 h L ∞ 1 + e n B (h) h d e n B (h) h f 0 L ∞ (50)
with a constant independent of h, f 0 , u and n. Moreover, if the exact and approximate backward flows satisfy uniform Lipschitz estimates

| Bn ex | 1 , sup k∈Z d | Bn h,k | 1,Σ n h,k ≤ C, (51)
then the transport error is also controlled by the forward flow error,

( T n h -T n ex )f 0 h L ∞ 1 + e n F (h) h d e n F (h) h f 0 L ∞ . (52)
Proof For x ∈ R d , we now denote

K n ex (x) := {k ∈ Z d : x ∈ F n ex (Σ 0 h,k )} and K n h (x) := {k ∈ Z d : x ∈ F n h,k (Σ 0 h,k )}.
The cardinality of K n ex (x) is readily bounded by the overlapping of the initial supports Σ 0 h,k : from ( 16) we find indeed

# K n ex (x) ≤ (2ρ 0 ) d . Moreover, for k ∈ K n h (x) we write hk -Bn ex (x) ∞ ≤ hk -Bn h,k (x) ∞ + Bn h,k (x) -Bn ex (x) ∞ < hρ 0 + e n B (h)
and since K n h (x) is a subset of Z d , the above bound yields

# K n h (x) ≤ 2(ρ 0 + h -1 e n B (h)) d .
It follows that the overlapping constant (45) is bounded by

Θ n (h) ≤ sup x∈R d #(K n h (x)) + #(K n ex (x)) (1 + h -1 e n B (h)) d , (53)
so that Lemma 1 gives the first estimate (50). In order to derive an estimate based on the forward flow error we next write for

k ∈ K n ex (x) that Bn h,k (x) -Bn ex (x) ∞ ≤ Bn h,k F n ex ( Bn ex (x)) -Bn h,k F n h,k ( Bn ex (x)) ∞ ≤ | Bn h,k | 1,Σ n h,k F n ex ( Bn ex (x)) -F n h,k ( Bn ex (x)) ∞ ≤ | Bn h,k | 1,Σ n h,k e n F (h)
where we have used that Bn ex (x) ∈ Σ 0 h,k in the last two inequalities. Similarly, for k ∈ K n (x) we write

Bn h,k (x) -Bn ex (x) ∞ ≤ Bn h,k (x) -Bn ex F n h,k ( Bn h,k (x)) ∞ ≤ | Bn ex | 1 F n ex ( Bn h,k (x)) -F n h,k ( Bn h,k (x)) ∞ ≤ | Bn ex | 1 e n F (h)
where we have now used that Bn h,k (x) ∈ Σ 0 h,k . According to (46) we thus have

e n B (h) ≤ max | Bn ex | 1 , sup k∈Z d | Bn h,k | 1,Σ n h,k e n F (h)
which gives the desired estimate. ⊓ ⊔

Particle transport with polynomial shape transformations

In order to describe particle methods with shape transformation of arbitrary orders and establish error estimates, we now assume that the exact flow F n ex can be applied exactly as well as its derivatives. Thus, in the traditional method the particles keep their shape and are simply translated with

T n h,(0) ϕ 0 h,k (x) = ϕ 0 h,k ( Bn (0),k (x)) = ϕ h (x -F n ex (x 0 k )), (54)
which corresponds to approximating the exact backward flow with

Bn (0),k (x) := x -F n ex (x 0 k ) + x 0 k .
For point (Dirac) particles this operator coincides with ( 14) and is exact. For finitesize particles however, the method does not converge in general. Assume indeed that ϕ is the hat function, and consider the smooth 2d problem where f 0 = 1 and u(t, x) = (-x 2 , x 1 ) over the time interval [0, π 4 ]. Then any reasonable initialization will give w k = h 2 , hence, f 0 h (x) = 1, and clearly the exact final solution is f ( π 4 , x) = 1. Now, at the final time the particle centers will have rotated of T = π 4 , therefore every particle with |k 1 | + |k 2 | = 1 will be centered on (cos(θ + π 4 ), sin(θ + π 4 )) with θ ∈ π 2 N, and hence contributes to x = 0 with T n h,(0

) ϕ 0 h,k (0) = h -2 (1 -1 √ 2 )
2 , in addition to ϕ h,0 which does not move. Since the other particles do not contribute to x = 0, the final error satisfies

( T n h,(0) -T n ex )f 0 h L ∞ ≥ | T n h,(0) f 0 h (0) -1| = 2( √ 2 -1) 2 , regardless of h.
To improve the accuracy of the transport operator, the error estimates in Section 4.1 suggest to use higher-order approximations of the backward flow. Letting indeed

φ k (s) = φ k (s; x) := ( Bn ex -I)( F n ex (x 0 k ) + s(x -F n ex (x 0 k ))), ( 55
)
we see that the approximation Bn (0),k (x) ≈ Bn ex (x) corresponds to the lowest-order expansion φ k (0) ≈ φ k (1). We may then consider r-th degree expansions, r ≥ 1,

Bn (r),k (x) := x -F n ex (x 0 k ) + x 0 k + φ ′ k (0) + • • • + 1 r! φ (r) k (0) ≈ Bn ex (x). ( 56
)
Here we could have used the alternate φk (s) :

= (I -F n ex )(x 0 k + s( Bn ex (x) -x 0 k )) since φk (1) -φk (0) = φ k (1) -φ k (0)
, but we observe that the form of (55) gives

φ (r) k (s) = d l 1 •••l r =1 ∂ l 1 • • • ∂ l r ( Bn ex -I)( F n ex (x 0 k ) + s(x -F n ex (x 0 k ))) r i=1 (x -F n ex (x 0 k )) l i
(57) so that Bn (r),k is a polynomial mapping which coefficients involve derivatives of Bn ex at F n ex (x 0 k ), which can be written in terms of the derivatives of F n ex at x 0 k . Moreover, (57) allows to specify the accuracy of the Taylor expansions (56). Indeed for every

x in a localized domain ω ⊂ B ℓ ∞ ( F n ex (x 0 k ), hρ) with ρ > 0, we have Bn (r),k (x)-Bn ex (x) ∞ = 1 0 (1 -s) r r! φ (r+1) k (s) ds ∞ ≤ h r+1 ρ r+1 (r + 1)! | Bn ex | r+1, ω (58
) where ω denotes the convex hull of ω.

For r = 1, observing that J Bn ex ( F n ex (x 0 k )) = (J F n ex (x 0 k )) -1 we obtain Bn (1),k (x) = x 0 k + (J n k ) -1 (x -F n ex (x 0 k )) with J n k := J F n ex (x 0 k ),
so that the linearly-transformed particle (LTP) transport operator reads

T n h,(1) ϕ 0 h,k (x) := ϕ 0 h,k ( Bn (1),k (x)) = ϕ h (J n k ) -1 (x -F n ex (x 0 k )) . ( 59
)
As already noted, this corresponds to using for the k-th particle the exact transport operator associated with the linearized flow at

x 0 k , F n (1),k (x) := ( Bn (1),k ) -1 (x) = F n ex (x 0 k ) + J n k (x -x 0 k ).
We are thus in position to use Lemma 2.

Theorem 1 The LTP transport operator (59) satisfies

( T n h,(1) -T n ex )A h f 0 L ∞ hc n F 1 + hc n F d f 0 L ∞ (60) with c n F := | F n ex | 2 1 | Bn ex | 2
, and an unspecified constant depending only on p and d.

Proof Applying Lemma 2, we obtain

( T n h,(1) -T n ex )A h f 0 L ∞ 1 + e n B,(1) (h) h d e n B,(1) (h) h f 0 L ∞ (61)
where we have set

e n B,(1) (h) := sup k∈Z d Bn (1),k -Bn ex L ∞ (Σ n h,k ) (62)
and 16) one easily sees that the transported particles are supported on

Σ n h,k := F n (1),k (Σ 0 h,k ) ∪ F n ex (Σ 0 h,k ). Next from (
F n (1),k (Σ 0 h,k ) = F n ex (x 0 k ) + J n k (B ℓ ∞ (0, hρ 0 )) ⊂ B ℓ ∞ ( F n ex (x 0 k ), hρ 0 J n k ∞ )
Moreover, the supports of the exactly transported particles satisfy

F n ex (Σ 0 h,k ) = F n ex (B ℓ ∞ (x 0 k , hρ 0 )) ⊂ B ℓ ∞ ( F n ex (x 0 k ), hρ 0 | F n ex | 1 ). ( 63
) Using next J n k ∞ ≤ | F n ex | 1 we obtain Σ n h,k ⊂ B ℓ ∞ ( F n ex (x 0 k ), hρ 0 | F n ex | 1
) and estimate (58) yields

e n B,(1) (h) ≤ h 2 (ρ 0 | F n ex | 1 ) 2 2 | Bn ex | 2 h 2 c n F ( 64
)
with a constant depending only on p. This completes the proof.

⊓ ⊔

Remark 2 From (53) we also see that the particles transported with the LTP operator (59) have a bounded overlapping constant sup

x∈R d # {k ∈ Z d : T n h,(1) ϕ 0 h,k (x) = 0} (1 + hc n F ) d (65) with again c n F := | F n ex | 2 1 | Bn ex | 2
and a constant depending only on p and d.

When r > 1, we have seen that the support of ϕ 0 h,k ( Bn (r),k (•)) has no reason to be contained in a ball of radius O(h), and for this reason the transported particles need to be restricted to a priori domains as in (43). For proving convergence rates one may define these domains by transporting small extensions of the initial supports with the linearized forward flow. Specifically, we consider (using now a tilde to distinguish this set from the one in (62))

Σn h,k := F n (1),k (B ℓ ∞ (x 0 k , hρ n h,k )) with ρn h,k := ρ 0 + 1 h Bn (1),k -Bn ex L ∞ ( F n ex (Σ 0 h,k )) . ( 66
)
Theorem 2 For r ≥ 1, the r-th order transport operator defined by

T n h,(r) ϕ 0 h,k (x) := 1 Σn h,k (x)ϕ 0 h,k ( Bn (r),k (x)) (67) satisfies ( T n h,(r) -T n ex )A h f 0 L ∞ h r cn F,(r) (1 + hc n F,( 1
)
) d f 0 L ∞ (68)
with an unspecified constant independent of h, f 0 , u and n, and where for r ≥ 1 we have set cn

F,(r) = (ρ n | F n ex | 1 ) r+1 | Bn ex | r+1 and ρn = ρ 0 (1 + hρ 0 2 | F n ex | 2 1 | Bn ex | 2 ).
Proof We first check that the support (66) contains F n ex (Σ 0 h,k ). To do so we take

x = F n ex (x) with x ∈ Σ 0 h,k and write Bn (1),k (x) -x 0 k ∞ ≤ ( Bn (1),k -Bn ex )(x) ∞ + x -x 0 k ∞ ≤ hρ n h,k .
This shows that x ∈ Σn h,k , hence Assumption (47) is fulfilled indeed and Lemma 1 applies: For some constant independent of h, f 0 , u and n, we have

( T n h,(r) -T n ex )A h f 0 L ∞ ẽn B,(r) (h) Θn (h) h f 0 L ∞ (69)
with an overlapping constant defined similarly as in (45), i.e., Θn (h) := sup

x∈R d # {k ∈ Z d : x ∈ Σn h,k } (70)
and a backward flow error defined similarly as in (46), i.e., ẽn B,(r) (h) := sup

k∈Z d Bn (r),k -Bn ex L ∞ ( Σn h,k ) .
To further bound the overlapping constant we proceed similarly as in the proof of Lemma 2: given x and k such that x ∈ Σn h,k , we write

hk -Bn ex (x) ∞ ≤ x 0 k -Bn (1),k (x) ∞ + Bn (1),k (x) -Bn ex (x) ∞ < hρ n h,k + ẽn B,(1) (h) ≤ hρ 0 + 2ẽ n B, (1) (h)
, and using that k ∈ Z d we find

Θn (h) ≤ 2ρ 0 + 4 ẽn B,(1) (h) h d . (71)
It remains to estimate the flow errors. From (66) and ( 62)-(64) we first derive

ρn h,k ≤ ρ 0 + e n B, (1) (h)
h ≤ ρ 0 + h 2 (ρ 0 | F n ex | 1 ) 2 | Bn ex | 2 = ρn . (72)
We then observe that any

x ∈ Σn h,k reads x = F n (1),k (x) = F n ex (x 0 k ) + J n k (x -x 0 k ) for some x ∈ B ℓ ∞ (x 0 k , hρ n h,k ). Using (72) and J n k ∞ ≤ | F n ex | 1 this gives Σn h,k ⊂ B ℓ ∞ ( F n ex (x 0 k ), hρ n h,k J n k ∞ ) ⊂ B ℓ ∞ ( F n ex (x 0 k ), hρ n | F n ex | 1 ).
Thus we can apply (58) with ω = Σn h,k , and take the supremum over k: this yields

ẽn B,(r) (h) h r+1 (ρ n | F n ex | 1 ) r+1 | Bn ex | r+1 = h r+1 cn F,(r) (73)
which completes the proof. ⊓ ⊔ Remark 3 (heterogeneous particles) Since our analysis does not rely on a smoothing kernel argument it readily extends to "particles" ϕ 0 h,k that are not derived from a reference function ϕ. In particular, our method applies to the transport of continuous finite element functions defined on an unstructured mesh, and the same results hold true under the usual regularity and quasi-uniformity assumptions.

Remark 4 (extension to non-linear problems) A important extension to the present work is the case where the (generalized) velocity field u is not given but instead depends on the transported density f . A typical example is the 1d1v Vlasov-Poisson system where u(t, x, v) = (v, -E(t, x)) and ∂ x E(t, x) = -f (t, x, v) dv. Here if the particle centers are pushed forward with a standard leap-frog scheme

x n+1 k = x n k + ∆tv n+1/2 k , v n+1/2 k = v n-1/2 k -∆tE n (x n k ), the corresponding forward flow reads F n : (x, v) → (x + ∆tv -∆t 2 E n (x), v -∆tE n (x)
) and an incremental deformation matrix is updated according to [START_REF] Raviart | An analysis of particle methods[END_REF] with

D n+1 k := D n k (J n k ) -1 and J n k ≈ J F n (x n k , v n-1/2 k ) = 1 -∆t 2 (E n ) ′ (x n k ) ∆t -∆t(E n ) ′ (x n k ) 1 .
This case has been numerically tested in [START_REF] Pinto | Noiseless Vlasov-Poisson simulations with linearly transformed particles[END_REF] and a priori error estimates were derived in [START_REF] Pinto | Uniform convergence of a linearly transformed particle method for the Vlasov-Poisson system[END_REF].

Numerical experiments

In this section we test the proposed LTP and QTP schemes and compare them with the standard TSP and FSL methods described in Section 2, using either M ′ 4 or cubic B-spline particles as described in Section 3.1. To assess the robustness of the method with respect to the velocity field and the initial data we use several passive transport problems in 2d, see Table 1. For the velocity fields we consider

• the reversible "swirling" velocity field proposed by LeVeque [START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF],

u SW (t, x; T ) := cos πt T curl φ SW (x) with φ SW (x) := -sin 2 (πx 1 ) sin 2 (πx 2 ) π • another reversible velocity field emulating a Raylegh-Benard convection cell,

u RB (t, x; T ) := cos πt T curl φ RB (x) with φ RB (x) := x 1 -1 2 (x 1 -x 2 1 )(x 2 -x 2
2 ) ; • and finally a constant non-linear rotation field derived from Example 2 in [START_REF] Bokanowski | An adaptive sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman equations[END_REF],

u NLR (x) := α(x) 1 2 -x 2 x 1 -1 2 with α(x) := 1 - x -( 1 2 , 1 2 ) 2 0.4 3 + .
Here the form of u SW and u RB yields reversible problems: at t = T /2 the solutions reach a maximum stretching, and they revert to their initial value at t = T . As for the non-linear rotation field u NLR , it is associated with the exact backward flow

Bn ex (x) = 1 2 1 2 + cos(α(x)t n ) sin(α(x)t n ) -sin(α(x)t n ) cos(α(x)t n ) x 1 -1 2 x 2 -1 2 ,
and the exact solutions are given by f (t n , x) = f 0 ( Bn ex (x)). In addition to the above velocity fields we consider the following initial data:

• smooth humps of approximate radius 0.2 given by

f 0 hump (x; x) := 1 2 1 + erf 1 3 (11 -100 x -x 2 )
and centered on x = (0.5, 0.4) or (0.5, 0.7), depending on the cases ; • a cone of radius 0.15 centered on x = (0.5, 0.25), f 0 cone (x; x) := 1 -20 3 x -x 2 + • and finally for the non-linear rotation field u NLR we take an initial data corresponding to Example 2 from [START_REF] Bokanowski | An adaptive sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman equations[END_REF], i.e.,

f 0 (x) := x 2 -1 2
. By combining the above values we obtain the four test-cases defined in Table 1, and accurate solutions are shown in Figures 1234for the purpose of illustration. In Table 1 we also give the respective time steps ∆t used in the time integration of the particle trajectories. In every case indeed, the numerical flow F n is computed with a RK4 scheme, and the time steps have been taken small enough to have no significant effect on the final accuracy. It happens that in every case we have ∆t = T /100, but this is unintended.

u(t, x) f 0 (x) T ∆t SW-cone u SW (t, x; T ) f 0 cone (x; x) with x = (0.5, 0.25) 5 0.05 SW-hump u SW (t, x; T ) f 0 hump (x; x) with x = (0.5, 0.7) 5 0.05 RB-hump u RB (t, x; T ) f 0 hump (x; x) with x = (0.5, 0.4) 3 0.03 NLR u NLR (x) x 2 -1 2 50 0.5

Numerical convergence rates

To measure the convergence properties of the proposed methods, we plot in Figures 1-4 the relative L ∞ errors at t = T versus the average number of active particles. Results are shown for the four test-cases defined in Table 1, and for each case the LTP and QTP methods are compared with the standard TSP and FSL methods described in Section 2, using either M ′ 4 or cubic B-spline particles. For the reversible test-cases we always remap the particles at t = T /2 when the solutions are stretched most, in order to take into account the accuracy of the intermediate approximations in the final measurements.

For every method we also compare different runs obtained by varying their main parameter: with the TSP method we take different values of the exponent q for which the particle scale ε behaves like h q (smaller values of q corresponding to more particle overlapping, see Section 2.1), and with the remapped particle schemes we test different values of the remapping period ∆t r = N r ∆t. From these results we make the following observations.

• As expected from theory, the TSP runs (shown on the first rows) only converge for values of q < 1, which corresponds to an extended overlapping: the ratio ε/h must go to +∞ as h → 0. We also observe that the convergence is always slow when compared to the remapped particle methods, despite the fact that the M ′ 4 kernel satisfies a fourth-order moment condition, see Estimate [START_REF] Pinto | Uniform convergence of a linearly transformed particle method for the Vlasov-Poisson system[END_REF]. • In order to converge, the FSL method (second rows) must be run with very short remapping periods, making it somehow closer to an Eulerian method. When remapped every few time steps indeed, the FSL runs exhibit significantly faster convergence rates than the TSP method. However, for fixed values of ∆t r there is always a point where the convergence stops. This amounts to asking for more remappings when more particles are used, which in practice has an effect similar to imposing a CFL constraint. • The most striking results with the LTP and QTP runs is that they seem to completely suppress the loss of convergence observed with the FSL method. Moreover, the observed behavior is now radically different: not only does the convergence always hold, it is improved when the remapping period is increased up to a certain value that is significantly larger than the time step. Specifically, we observe that the accuracy of the LTP runs improves for remapping periods as large as 10 to 50 times the time step, depending on the test cases. And with the QTP method these ratios go up to values between 30 and 50. • Broadly speaking, the above hold true for both the M ′ 4 and the B 3 particles, the main difference being that the former perform better for small remapping periods. Given the fact that they are remapped with lower-order but also lower dissipation than the B 3 particles, this behavior is rather expected.

The remaining sections are devoted to further investigating the influence of the remapping period ∆t r on the accuracy.

Influence of the remapping period

In Figure 5 we plot the final errors obtained with different runs using increasing values for the remapping period ∆t r . This leads us to the following observations.

• With the FSL scheme the particles must be remapped each few time steps, as the accuracy of the method quickly deteriorates for increasing remapping periods. This confirms the previous observations. • By transforming the particles either linearly or quadratically, we obtain a twofold benefit: first, the accuracy is always improved, and sometimes significantly. Second, the optimal accuracy now corresponds to some trade off between small remapping periods where the remappings error dominate the transport errors, and large ones where the opposite occurs. • By transforming the particles quadratically instead of linearly, we observe some further improvements in the accuracy, but the major benefit seems to be a significant gain in robustness with respect to the remapping period. In the NLR test-case for instance, the optimal remapping period for the QTP scheme is about five time larger than for the LTP scheme, be it with the moderate runs shown in Figure 5, or with finer ones using 512 × 512 particles, not shown here.

• Finally, we verify that neither the LTP nor the QTP schemes are subject to a CFL-like condition. Indeed, if that was the case the optimal remapping periods should decrease as h → 0, but we observe the opposite behavior.

Numerical study of the dynamic remapping criterion

In Figure 6 we plot the results obtained with the dynamic criterion (39) and compare them with those obtained with a static strategy. To do so, for each case we plot both error curves using as x-axis an average remapping period : for the static runs it is the constant remapping period (hence the curves correspond to those already shown in Figure 5) and for the dynamic runs it is defined as T /(1+R) where R is the number of dynamic remappings (initialization included).

Finally, the different points in the dynamic runs correspond to different values of the constant in (39). For the LTP runs the red (resp. blue) points correspond to values larger (resp. smaller) than C remap = 1. For the QTP runs the threshold value is C remap = 5.

From the results shown in Figure 6 we may draw a positive conclusion: indeed in every case but one, our dynamic remapping strategy achieves the same level of accuracy as the best static run (in some cases it is even more accurate), while using about the same number of remappings. And in the NLR case where the dynamic strategy fails to reach that optimal accuracy one could argue that the solutions are very smooth, hence by only considering the first order error estimate (41) we obtain a remapping error indicator that is overly pessimistic. One would expect better results with a numerical indicator based on a second order estimate.

Conclusion

In this work we have proposed and analyzed discrete implementations of a new class of particle methods where the shape functions are transformed to better approximate the transport along the characteristic flow.

To reduce the smoothing involved in the most basic particle methods which typically involve convolutions by particles with extended overlapping or very frequent remappings, we have first described a linearly-transformed particle (LTP) scheme and we have shown that it achieves good uniform accuracy with bounded particle overlapping and low remapping frequencies compared to fixed-shape particle methods.

Because the linearization process underlying the LTP method amounts to locally discard the second-order derivatives of the characteristic flow, which one may see as a smoothing by itself, we also proposed a quadratically-transformed particle (QTP) scheme where this effect is further reduced. By a numerical study involving several 2d test problems we have shown that this approach allowed to further reduce the optimal remapping frequencies by a significative amount.

On a theoretical level, we have extended these methods to polynomial shape transformations of arbitrary degree r, and for problems with smooth characteristic flow we have established their L ∞ convergence with order r. As our estimates require no particle remappings, they support the observed robustness of the methods with respect to the remapping frequency. Specifically, we have observed that with the LTP scheme the optimal remapping periods could be as large as about 10 to 50 times what they are for standard (fixed-shape) particle methods, and with the QTP scheme this ratio often reaches values beyond 30.

Finally, we have proposed a dynamic criterion to only remap the particles when needed. Our strategy is supported by the error analysis, and it is validated by the numerical tests.

An important feature of the proposed schemes is that they only involve pointwise evaluations of the forward flow: since virtually every particle code contains routines that compute accurate approximations of this flow to push the particles forward in time, we believe that our schemes should be simple to implement within existing codes, with reasonable additional programming cost. This property has been tested in a recent work [START_REF] Pinto | Noiseless Vlasov-Poisson simulations with linearly transformed particles[END_REF] where an LTPIC scheme has been studied for a Vlasov-Poisson plasma. Fig. 5 (Color) Relative L ∞ errors at t = T vs. remapping period ∆tr for the different test cases solved with remapped particle methods of order 0 (FSL), 1 (LTP) and 2 (QTP). Here the particles are initialized and remapped on a cartesian grid with h = 2 -8 , which correspond to a maximum of 256 × 256 particles (in some cases only a fraction of these particles are activated). Left panels show results obtained with M ′ 4 particles, whereas right panels show results obtained with cubic B-spline particles. Qualitatively similar curves have been obtained with runs using finer grids of 512 × 512 particles, not shown here. Fig. 6 (Color) Performance of the dynamic remapping strategy applied to the LTP and QTP method. Here the relative L ∞ errors at t = T are plotted vs. the average remapping periods, for the four test cases defined in Table 1 (every run is obtained with 256 × 256 particles). In the static runs the abscissa represents the constant remapping period ∆tr as in Figure 5.

In the dynamic runs it corresponds to the ratio T /R where R is the number of remappings (initialization included) selected by the criterion (39) with the error indicators (40), (42). Finally the different points in the dynamic runs correspond to different values of the constant in (39). For the LTP runs the red (resp. blue) points correspond to values larger (resp. smaller) than Cremap = 1. For the QTP runs the threshold value is Cremap = 5.

5 Fig. 1 ( 5 Fig. 2 ( 5 Fig. 3 ( 10 Fig. 4 (

 515253104 Fig.1(Color) Convergence curves (relative L ∞ errors at t = T vs. average number of active particles) for the reversible test case SW-cone defined in Table1, solved with the different methods (see text for details). The first row shows the profile of the exact solution: the initial (and final) density f 0 = f (T ) is on the left, whereas the intermediate solution f (T /2) (with maximum stretching) is on the right.

Table 1

 1 Definition of the benchmark test-cases name

Acknowledgements The author thanks Eric Sonnendrücker, Albert Cohen, Jean-Marie Mirebeau and Jean Roux for valuable discussions during the early stages of this research.


François Sarrazin

Philippe Pouliguen

Ala Sharaiha

Patrick Potier

Janic Chauveau

Window increasing technique to discriminate mathematical and physical resonant poles extracted from antenna response

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The window increasing technique to discriminate mathematical and physical resonant poles extracted from antenna response F. Sarrazin, P. Pouliguen, A. Sharaiha, P. Potier and J. Chauveau This letter presents a new approach called window increasing technique (WIT) to discriminate mathematical and physical poles extracted from a noisy antenna response. The principle of the WIT is to apply a pole extraction method on several windows of the response and then to observe the stability of the extracted poles. In order to compare the WIT to the classical window moving technique (WMT), we apply these two techniques on the electric far field backscattered by a dipole antenna. We show that, in presence of noise, the WIT allows finding more physical poles with a good accuracy than the WMT.

Since its introduction by Baum [1], the singularity expansion method (SEM) has been widely used for antenna characterisation [2/3]. This method allows modelling the late time response of an antenna with only a few sets of parameters: poles and residues. The main advantage of the SEM is that resonant poles, also called complex natural resonance (CNR), depend only on antenna characteristics and are independent of the direction of the incoming wave, the excitation waveform and its polarization [1]. Therefore, poles allow representing an antenna in a compact unique way. There are several extraction methods to obtain CNR from antenna response but the most commonly used is the matrix pencil method (MPM) [4]. In practice, the number of poles contained in a response is unknown; that is why the number of poles to be extracted is usually overestimated by the MPM. Due to this overestimation, some mathematical poles are extracted in addition to the physical ones. In order to use the CNR to characterize an antenna, one needs to discriminate the mathematical poles from the physical ones. A classical way to discriminate poles is the window moving technique (WMT). However, the WMT is limited when applied on noisy response [5]. In this letter, we suggest a new approach called the window increasing technique (WIT). In order to compare the WMT and the WIT, these two techniques are applied on the noisy electric far field backscattered by a dipole antenna.

We consider a dipole antenna of length 34 mm and diameter 0.05 mm, so its ratio / is equal to 680. The impedance of its lumped port is a matched load of 73 =. This antenna, simulated using CST Microwave Studio, is excited by a plane wave in the boresight direction and the backscattered electric far field is measured using a probe at a distance of 2 meters. The late time response backscattered by the dipole antenna is shown in Fig. 1. For the noisy case, a white Gaussian noise (WGN) is added to the noiseless response to obtain a signal to noise ratio (SNR) of 10 dB.

The MPM is applied to the noiseless response and the resonant poles are shown in the complex plane in Fig. 2 with a | | | | ⁄ weighting where is the damping coefficient of the pole, i.e. its real part and is its residue. It means that the more the marker is big, the more the pole's contribution is important. Resonant frequencies of poles correspond to natural frequencies of the dipole antenna at /2, 3 /2, 5 /2 et 7 /2 where λ is the free space wavelength.

| | | | ⁄

The principle of the WMT is to apply the MPM on a windowed response. Then, the window is shifted of a small time step and the MPM is applied again. The assumption is that, depending on the window, the position of the mathematical poles will change from window to window whereas the physical poles will remain essentially unchanged. The WMT is applied on the noisy late time dipole response shown in Fig. 1. The window length is 140 samples (0.29 ns) and the window shift is 10 samples (0.019 ns). Results are presented in terms of resonant frequencies and damping coefficients in Fig. 3. We can see that three resonant frequencies are stable around 4, 13 and 21 GHz. These frequencies are in good agreement with frequencies of the three first poles extracted from the noiseless response. Moreover, the WMT confirms the poles' weighting. Indeed, the more a pole's weight is important, the more this pole is accurately extracted late from the response. However, damping coefficients are not as stable as resonant frequencies. The two first damping coefficients can be evaluated around 3. 10 and 4. 10 Neper/s but the third one cannot be defined. Indeed, its variation is too important according to the window. Using this approach, only two physical poles of the dipole antenna can be determined with a good accuracy.

!"

Resonant frequencies

Since results of the WMT seem to be stable for the first windows, we suggest keeping the beginning of the antenna response in all windows considered. Indeed, these first samples contain data with the highest SNR. This novel approach, that we call the WIT, consists of applying the MPM to a windowed response beginning from its first samples, and then to increase the window's length and to apply the MPM again. It means that the beginning time of the window is unchanged whereas the end time is increased until including the final sample. Results of the WIT applied on the noisy dipole antenna response are presented in terms of resonant frequencies and damping coefficients in Fig. 4.

!

Resonant frequencies Damping coefficients

Three resonant frequencies around 4, 13 and 21 GHz are perfectly stable for all windows and the fourth one is stable around 31 GHz except for some windows. All these frequencies are in good agreement with ones extracted from noiseless response. Moreover, damping coefficients of the three first poles are also very stable according to the window and correspond to those extracted in the noiseless case.

#

In this letter, the novel WIT is proposed in order to discriminate the physical poles of an antenna response. This technique is compared to the classical WMT. Although the WMT could be of interest, the WIT appears to be a more powerful technique that allows extracting stable damping coefficient, especially in presence of noise.

$

This work was supported by the Direction Générale de l'Armement (DGA), from the French ministry of defence.

  F. Sarrazin and A. Sharaiha (Institute of Electronics and Telecommunications of Rennes, University of Rennes 1, France) Hua, Y. and Sarkar, T.K.: 'Generalized pencil/of/function method for extracting poles of an EM system from its transient response', %%%

								E/mail: francois.sarrazin@univ/rennes1.fr
								P. Pouliguen, P. Potier and J. Chauveau (Direction Générale de
								l'Armement, France)
								1 Baum, C.E.: 'On the singularity expansion method for the solution of
								electromagnetic interaction problems,' EMP Interaction Note 8, Air
								Force Weapons Laboratory, Kirkland AFB, New Mexico, 1971
								2 Sarrazin, F., Chauveau, J. Pouliguen, P. Potier, P. and Sharaiha, A.:
								'Accuracy of singularity expansion method in time and frequency
								domains to characterize antennas in presence of noise', for coming in
			35					%%%	$	&
	Resonant Frequency (GHz)	5 10 15 20 25 30					3 Marchais, C., Uguen, B., Sharaiha, A., Le Ray, G. and Le Coq, L.: 'Compact characterisation of ultrawideband antenna responses from frequency measurements '. % " $ & , 2011, , (6), pp. 671-675 4 $ & , 1989, , (2), pp. 229-234
			0.2 0	0.4	0.6 Window end time (ns) 0.8	1	1.2	5 Manteghi, M. and Rezaiesarlak, R.: 'Short time matrix pencil for chipless RFID detection applications', %%% $
								&	, 2013, , (5), pp. 2801-2806
	Neper/s)		-2 0
	Damping Coefficient (10 9	0.2 -10 -8 -6 -4	0.4	0.6	0.8	1	1.2
					Window end time (ns)