Abstract-The graph edit distance (GED) measures the amount of distortion needed to transform a graph into another graph. Such a distance, developed in the context of error-tolerant graph matching, is one of the most flexible tool used in structural pattern recognition. However, the computation of the exact GED is NP-complete. Hence several suboptimal solutions, such as the ones based on bipartite assignments with edition, have been proposed. In this paper we propose a binary quadratic programming problem whose global minimum corresponds to the exact GED. This problem is interpreted as a quadratic assignment problem (QAP) where some constraints on solutions have been relaxed. This allows to adapt the integer projected fixed point algorithm, initially designed for the QAP, to efficiently compute an approximate GED by finding an interesting local minimum. Experiments show that our method remains quite close to the exact GED for datasets composed of small graphs, while keeping low execution times on datasets composed of larger graphs.

I. INTRODUCTION

Computing an efficient similarity or dissimilarity measure between graphs is a major problem in structural pattern recognition. The graph edit distance (GED) [START_REF] Tsai | Error-correcting isomorphisms of attributed relational graphs for pattern analysis[END_REF]- [START_REF] Bunke | Error correcting graph matching: on the influence of the underlying cost function[END_REF], developed in the context of error-tolerant graph matching, provides such a measure. It may be understood as the minimal amount of distortion required to transform a graph into another, by a sequence of edit operations applied on nodes and on edges, restricted here to substitutions, insertions and removals. Such a sequence is called an edit path. Depending on the nature of the graphs and on the context, all paths may not have a same importance, in particular when nodes or edges are labeled or attributed. So each possible edit operation e is penalized by a non-negative cost c(e), and the integration of these costs over an edit path γ defines the length A(γ) of this path. An edit path having a minimal length, among all edit paths Γ(G 1 , G 2 ) transforming a graph G 1 into a graph G 2 , defines the GED from

G 1 to G 2 : GED(G 1 , G 2 ) = min γ∈Γ(G1,G2) A(γ) def. = e⊂γ c(e) (1)
Since computing the GED is NP-complete, it is restricted to rather small graphs. So several approaches have been proposed to approximate the GED efficiently and to process larger graphs. In this paper, graphs are assumed to be simple (no loop nor multiple edge), and each element of the two graphs can be edited only once (no composition of edit operations). Under these conditions, the GED can be expressed as a binary quadratic program (QP), which is the subject of this paper. Contrary to the graph matching problem, which can be rewritten as a quadratic assignment problem (QAP) [START_REF] Burkard | Assignment Problems[END_REF], [START_REF] Leordeanu | An integer projected fixed point method for graph matching and map inference[END_REF], an error-tolerant graph matching can also remove and insert elements. Managing these edit operations is the main difficulty of the quadratic framework. To this, a binary QP is obtained in [START_REF] Justice | A binary linear programming formulation of the graph edit distance[END_REF] by considering an edit grid. It is transformed into a binary linear program, which can be solved exactly with classical solvers, but graphs are restricted to be undirected, edges to be unlabeled, and it is also restricted to small graphs since there are (n + m) 2 binary variables, with n = |V 1 | and m = |V 2 |. In [START_REF] Neuhaus | A quadratic programming approach to the graph edit distance problem[END_REF] a QP with nm real variables is defined to optimize the cost restricted to node substitutions and associated edge operations. This QP, solved by an interior point algorithm, provides a weight for each possible node substitution. From this fuzzy assignment between nodes, an approximate minimal-cost edit path is then constructed using an heuristic method. Other approximate GED are based on a similar two-step strategy.

In particular, the bipartite graph edit distance (bGED) [START_REF] Riesen | Bipartite graph matching for computing the edit distance of graphs[END_REF], [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF] approximates the GED by replacing the error-tolerant matching between two graphs by an error-tolerant bipartite matching of their nodes, from which an edit path is constructed. Such an error-tolerant set matching can be rewritten as a square linear sum assignment problem (LSAP), with (n + m) 2 binary variables to take into account substitutions, insertions and removals of elements. To approximate the GED, these operations are penalized by the cost of substituting, removing or inserting substructures attached to the nodes [START_REF] Riesen | Bipartite graph matching for computing the edit distance of graphs[END_REF]- [START_REF] Riesen | Estimating graph edit distance using lower and upper bounds of bipartite approximations[END_REF]. Since the LSAP can be solved in polynomial time, for instance in O((n + m) 3 ) with the well-known Hungarian algorithm [START_REF] Burkard | Assignment Problems[END_REF], [START_REF] Lawler | Combinatorial Optimization: Networks and Matroids[END_REF], the bGED provides an efficient approximation of the GED. This linear framework has been extended to a quadratic one [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF], where the GED is expressed as a global solution of a QAP with (n + m) 2 binary variables. An interesting local minimal solution is obtained by relaxing the constraints on the solutions, and by using the integer projected fixed point algorithm (IPFP) [START_REF] Leordeanu | An integer projected fixed point method for graph matching and map inference[END_REF]. This algorithm iterates: an instance of the LSAP (projection of local continuous solution in the discrete domain), and a line search in the continuous domain to get the next local continuous solution.

While this last quadratic framework improves the approximation of the exact GED, the size of the initial data has been artificially augmented to manage removal and insertion operations more easily. This is also the case in the linear framework above, and in the one proposed by [START_REF] Justice | A binary linear programming formulation of the graph edit distance[END_REF]. Recently, the errortolerant set matching involved in the computation of bGEDs has been formulated as an extension of the LSAP with only (n + 1)(m + 1) binary variables [START_REF] Bougleux | Linear sum assignment with edition[END_REF]. Since a solution can be computed in O(min{n, m} 2 max{n, m}) time complexity by an adapted Hungarian algorithm, the computation of the bGED is improved, in particular for large graphs.

Based on this last formulation (Sec. II) [START_REF] Bougleux | Linear sum assignment with edition[END_REF], we propose in this paper a new quadratic expression of the GED with only (n + 1)(m + 1) binary variables (Sec. III). This problem is interpreted as a quadratic assignment problem (QAP) where constraints on solutions have been relaxed for removal and insertions operations. This allows to adapt the IPFP algorithm used in [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF] to compute efficiently an approximate GED. Experiments (Sec. IV) show that our method remains quite close to the exact GED for datasets composed of small graphs, as previously obtained in [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF]. Moreover, contrary to [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF], execution times on datasets composed of larger graphs are similar to the ones obtained by bGEDs.

II. LINEAR SUM ASSIGNMENT PROBLEM WITH EDITION

This section describes the problem of minimal-cost errortolerant set matching, formalized as an extension of the LSAP, as proposed in [START_REF] Bougleux | Linear sum assignment with edition[END_REF].

A. Assignment with edition

A transformation of a set V 1 into a set V 2 can be performed by applying edit operations as follows: Each element i ∈ V 1 is either substituted by an unique element j ∈ V 2 (denoted by i → j), or removed from V 1 (i → ). Then each remaining element j ∈ V 2 , not previously used for a substitution, is inserted (denoted by → j).

Consider the two sets

V 1 = V 1 ∪ { } and V 2 = V 2 ∪ { }.

Edit operations (substitutions, removals and insertions) defining the transformation can be represented by a mapping

ϕ : V 1 → P(V 2 ) satisfying:    ∀i ∈ V 1 , |ϕ(i)| = 1 ∀j ∈ V 2 , ϕ -1 [j] = 1 ∈ ϕ( ) (2)
where P(•) is the powerset, and ϕ -1 [j]

def.

= ϕ -1 [{j}] denotes the pre-image of any singleton {j} ∈ P(V 2 ) by ϕ.

Following [START_REF] Bougleux | Linear sum assignment with edition[END_REF], such a mapping is called an assignment with edition, or -assignment. It may be understood as a bijection on which the bijectivity constraint has been relaxed on the element . Indeed, we implicitly have

1 ≤ |ϕ( )| ≤ m + 1 and 1 ≤ |ϕ -1 [ ]| ≤ n + 1, with n = |V 1 | and m = |V 2 |. Let A (V 1 , V 2 ) be the set of -assignments from V 1 to V 2 .
To simplify the forthcoming expressions, consider w.l.o.g. that

V 1 = {1, . . . , n}, V 2 = {1, . . . , m}, V 1 = V 1 ∪ {n + 1} and V 2 = V 2 ∪ {m + 1}. Any -assignment ϕ ∈ A (V 1 , V 2 ) can be equivalently represented by a (n + 1) × (m + 1) binary matrix X = V 2         x 1,1 . . . x 1,m x 1,m+1 V 1 . . . . . . . . . . . . x n,1 . . . x n,m x n+1,m+1 x n+1,1 . . . x n+1,m 1 (3) such that        x i,j = δ ϕ(i)={j} , ∀(i, j) ∈ V 1 × V 2 (sub.) x i,m+1 = δ ϕ(i)={ } , ∀i ∈ V 1 (rem.) x n+1,j = δ ϕ -1 [j]={ } , ∀j ∈ V 2 (ins.) x n+1,m+1 = 1 (4)
Due to the constraints on ϕ, the matrix X has a 1 on each of its n first rows and a 1 on each of its m first columns:

   m+1 j=1 x i,j = 1, ∀i = 1, . . . , n n+1 i=1 x i,j = 1, ∀j = 1, . . . , m x n+1,m+1 = 1 (5)
Reciproqualy, any matrix X ∈ {0, 1} (n+1)×(m+1) satisfying Eq. 5 represents an -assignment.

B. Problem formulation

The definition of an -assignment does not rely on the nature of the underlying data. To select a relevant -assignment in A (V 1 , V 2 ), each edit operation o is penalized by a nonnegative cost c(o). Then the cost of an -assignment ϕ, is obtained by summing the costs of its edit operations:

L(ϕ) = i∈V 1 j∈ϕ(i) c(i → j) = i∈V1 ϕ(i)={j} c(i → j) substitutions + i∈V1 ϕ(i)={ } c(i → ) removals + j∈V2 ϕ -1 [j]={ } c( → j) insertions (6)
Note that the mapping → has a zero cost, since it does not represent any edit operation.

The cost of all edit operations, induced by

A (V 1 , V 2 ), can be collected into a matrix C ∈ [0, +∞) (n+1)×(m+1) such that c i,j = c(i → j) for all (i, j) ∈ V 1 × V 2 , with c n+1,m+1 = 0.
Consider an -assignment ϕ and its corresponding binary matrix X (Eq. 4). Then the cost of ϕ (Eq. 6) is equal to:

L(X) = n+1 i=1 m+1 j=1 c i,j x i,j = vec(C) T vec(X) (7)
where vec(M) is the vectorization of the p × q matrix M obtained by concatenating its rows, i.e. a vector of size pq.

A relevant -assignment is then defined as one having a minimal cost among all -assignments transforming V 1 into V 2 , i.e. satisfying argmin

ϕ∈A (V1,V2) L(ϕ) (8)
or equivalently satisfying the following binary linear program:

argmin x c T x | Lx = 1 n+m , x ∈ {0, 1} p , x p = 1 (9)
where p = (n + 1)(m + 1), c = vec(C) and the linear system Lx = 1 is the matrix version of the constraints defining an -assignment (Eq. 2 or Eq. 5). The assignment matrix L ∈ {0, 1} (n+m)×(n+1)(m+1) is defined by:

∀(i, j), l k,(i,j) = δ k=i , ∀k = 1, . . . , n l n+k,(i,j) = δ k=j , ∀k = 1, . . . , m (10)
Since L is totally unimodular, from standard tools in linear programming, the problem defined by Eq. 9 (or Eq. 8) has a binary solution. This problem is close to the classical linear sum assignment problem (LSAP). Compared to the LSAP, elements can also be removed or inserted, i.e. the constraints in Eq. 9 are relaxed for one element ( ) in each of the two sets. For this reason, the problem defined by Eq. 9 is called the linear sum assignment problem with edition (LSAPE) [START_REF] Bougleux | Linear sum assignment with edition[END_REF]. The LSAPE can be efficiently solved in polynomial time complexity by adapting algorithms that initially solve the LSAP. An adaptation of the Hungarian algorithm [START_REF] Burkard | Assignment Problems[END_REF], [START_REF] Lawler | Combinatorial Optimization: Networks and Matroids[END_REF] is presented in [START_REF] Bougleux | Linear sum assignment with edition[END_REF]. It computes a solution to the LSAPE in O(min{n, m} 2 max{n, m}) time complexity and in O(nm) space complexity. This improves complexities obtained by the classical approach used to transform a set into another with edit operations [START_REF] Riesen | Bipartite graph matching for computing the edit distance of graphs[END_REF], [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF]. In this last approach, the problem is formalized as a LSAP with a specific (n + m) × (m + n) cost matrix. It can be solved by the classical Hungarian algorithm in O((n + m) 3 ) time complexity and in O((n + m) 2 ) space complexity.

III. A QUADRATIC EXPRESSION OF THE GED

In this section we extend the linear sum assignment problem with edition to a quadratic one. This simplifies the quadratic expression of the GED proposed in [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF].

A. Simultaneous assignment with edition

Let G 1 = (V 1 , E 1 ) and G 2 = (V 2 , E 2 ) be two simple graphs, directed or not. A transformation of G 1 into G 2 ,
with edit operations applied on nodes and on edges, can be fully represented by an -assignment ϕ :

V 1 → P(V 2 ) transforming V 1 into V 2 .
Then the set E 1 is transformed into the set E 2 as follows:

1) Any edge (i, j) ∈ E 1 is mapped by ϕ to a pair (k, l) ∈ V 2 × V 2 with ϕ(i) = {k} and ϕ(j) = {l}. The edge

(i, j) is either a) substituted by (k, l) if (k, l) ∈ E 2 , or b) removed if (k, l) ∈ E 2 , in particular if i or j is removed by ϕ (k = or l = ). 2) Any edge (k, l) ∈ E 2 is inserted into E 1 if it has not
been used by a substitution in Step 1, i.e. if there is a pair

(i, j) ∈ V 1 × V 1 such that k ∈ ϕ(i), l ∈ ϕ(j) and (i, j) ∈ E 1 .
In particular, this occurs when k or l is inserted by ϕ into V 1 (i = or j = ). The sequence of edit operations transforming E 1 into E 2 is thus represented by a simultaneous -assignment of the nodes. Its cost is defined by summing the cost of its edit operations:

Q(ϕ) = (i,j)∈E1 with {k}=ϕ(i),{l}=ϕ(j),(k,l)∈E2 c((i, j) → (k, l)) (sub.
)

+ (i,j)∈E1 with {k}=ϕ(i),{l}=ϕ(j),(k,l) ∈E2 c((i, j) → ) (rem.) + (k,l)∈E2 with {i}=ϕ -1 [k],{j}=ϕ -1 [l],(i,j) ∈E1 c( → (k, l)) (ins.) (11)
i, j ∈ V 1 TO k, l ∈ V 2 RESPECTIVELY (i, j) (k, l) edit operation cost d ik,jl ∈ E 1 ∈ E 2 substitution of (i, j) by (k, l) c((i, j) → (k, l)) ∈ E 1 ∈ E 2 removal of (i, j) c((i, j) → ) ∈ E 1 ∈ E 2 insertion of (k, l) into E 1 c( → (k, l)) ∈ E 1 ∈ E 2 do nothing 0
More generally, the cost of simultaneously assigning i ∈ V 1 to k ∈ V 2 , and j ∈ V 1 to l ∈ V 2 , can be defined by:

d ik,jl = c ((i, j) → (k, l)) δ (i,j)∈E1 δ (k,l)∈E2 (sub.) + c((i, j) → )δ (i,j)∈E1 (1 -δ (k,l)∈E2 ) (rem.) + c( → (k, l))(1 -δ (i,j)∈E1 )δ (k,l)∈E2 (ins.) (12
) where δ r = 1 if the relation r is true or δ r = 0 else. The edit operation associated to the simultaneous node assignment is reported in Table I. Observe that d ik,jl = 0 when both (i, j) and (k, l) are not edges of G 1 and G 2 respectively. Indeed in this case the simultaneous assignment does not correspond to an edit operation applied on edges. Then the cost Q(ϕ) defined by Eq. 11 can be written more compactly as:

Q(ϕ) = g i∈V 1 k∈ϕ(i) j∈V 1 l∈ϕ(j) d ik,jl with g = 1 2
if G 1 and G 2 are undirected 1 else [START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF] When graphs are both undirected, the cost d ik,jl is symmetric (d ik,jl = d jl,ik ). So a factor 1 2 is introduced to count the cost of each edit operation only once.

B. Quadratic assignment problem with edition and GED

An optimal transformation of G 1 into G 2 consists in computing an -assignment having a minimal cost that penalizes operations both on nodes and edges, among all -assignments from V 1 onto V 2 . In other terms, the cost of this optimal transformation is equal to the GED from G 1 to G 2 [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF]:

GED(G 1 , G 2 ) = min ϕ∈A (V1,V2) A(ϕ) def. = Q(ϕ) + L(ϕ) (14)
This optimization problem is close to the quadratic assignment problem (QAP), see [START_REF] Burkard | Assignment Problems[END_REF] for more details on QAP.

To see this, consider an -assignment ϕ ∈ A (V 1 , V 2 ) and the associated binary vector x ∈ {0, 1} p , with x p = 1 and p = (n + 1)(m + 1). Nodes i and j of V 1 are simultaneously assigned by ϕ to k and l of V 2 , respectively, iff x ik = x jl = 1.

The quadratic cost Q(ϕ) can then be rewritten as

Q(x) = g n+1 i=1 m+1 k=1 n+1 j=1 m+1 l=1 d ik,jl x ik x jl = g x T Dx ( 15
)
where the matrix D ∈ [0, +∞) p×p , of general term d ik,jl (Eq. 12), encodes the cost of all possible simultaneous assignments. Compared to Eq. 13, the sum in Eq. 15 is also carried on simultaneous assignments which are not induced by the -assignment ϕ, i.e. x ik = 0 or x jl = 0. Since they do not contribute to the sum, we have Q(x) = Q(ϕ). Due to the constraints on -assignments, when i, j ∈ V 1 and k, l ∈ V 2 remark that several simultaneous assignments never occur, i.e. x ik x il = 0 with k = l or x ik x jk = 0 with i = j (a node cannot be assigned twice). In these cases, the value of d ik,jl is not well-defined by Eq. 12 and should be equal to 0. Nevertheless d ik,jl x ik x jl = 0 is always satisfied according to the previous remark, and so the value of the cost is not so relevant. Also, since graphs are simple, all diagonal elements are null (d ik,ik = 0). The cost of transforming G 1 into G 2 can then be written as

A(x) = g x T Dx + c T x ( 16
)
where c ∈ [0, +∞) p is the edit cost vector defined in Eq. 9.

Since x T Dx =1 2 x T (D + D T )x, we have

A(x) = 1 2 x T ∆x + c T x with ∆ = D if G 1 and G 2 are undirected D + D T else (17)
The cost matrix ∆ is always symmetric, even if at least one of the two graphs is directed. The linear part of the cost can also be included in the quadratic part, leading to

A(x) = x T ∆x, with ∆ = 1 2 ∆ + diag(c) (18)
Then the GED is rewritten as the following quadratic program:

GED(G 1 , G 2 ) = min x T ∆x | Lx = 1 n+m , x ∈ {0, 1} p (19)
We call it the quadratic assignment problem with edition (QAPE). The QAPE differs from the classical QAP by the constraint matrix L. A QAP finds a bijection or an injection inducing a simultaneous assignment, as for the quadratic expression of the GED proposed in [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF]. This expression describes a QAP involving assignments represented by vectors of size (n + m) 2 . This is simplified by Eq. 19 by involving vectors of size (n + 1)(m + 1) only.

C. Relaxation, approximation, algorithm

While exact solutions to the QAPE could be computed by adapting algorithms that initially solve the QAP [START_REF] Burkard | Assignment Problems[END_REF], this would not be competitive with bipartite GEDs in terms of computational time. So we propose to approximate the GED by adapting the integer projected fixed point (IPFP) algorithm [START_REF] Leordeanu | An integer projected fixed point method for graph matching and map inference[END_REF], designed to approximate the solution to the QAP, and it is also used in [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF] to approximate the GED. As several other algorithms, the IPFP algorithm tries to find a solution to the relaxed quadratic problem given here by

argmin x T ∆x | Lx = 1 n+m , x ∈ [0, 1] p (20)
where the space of solutions is now continuous.

Given an initial candidate solution x, the algorithm iterates the three following steps

b ← argmin (x T ∆)b | Lb = 1, b ∈ {0, 1} p , b p = 1 (21) α ← argmin α∈[0,1] A(x + α(b -x)) (22) x ← x + α (b -x) (23)
until convergence (when x becomes discrete) or after a given number of iterations. The first step (Eq. 21) consists in computing a linear approximation of A around the current solution x in the discrete domain. This can be realized by solving a LSAPE with x T ∆ as an edit cost matrix (Sec. II). Indeed, the firstorder Taylor expansion of A around x is given by

A(y) ≈ A(x) + (x T ∆)(y -x), y ∈ [0, 1] p (24)
so minimizing A(x) locally is approximatively equivalent to minimizing x T ∆y with x fixed. By standard tools in linear programming, this linear program has a same binary solution b as the problem defined by Eq. 21, here a LSAPE. It can be solved in O(min{n, m} 2 max{n, m}) time complexity and in O(nm) space complexity (Sec. II-B). Compared to the IPFP algorithm proposed in [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF], where Eq. 21 corresponds to a square LSAP, computational time is reduced at each iteration. This confirmed experimentally (Section IV). The second step consists in computing the local minimum of A, by a line search, between the current solution x and the candidate binary solution b included. As detailed in [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF] this can be solved analytically, here in O(mn). The vector realizing this minimum is the new candidate solution.

The iterative process generally converges to a local minimum of the relaxed problem defined by Eq. 20. When x is continuous, an additional projection (Eq. 21) is performed to get the final binary vector. The quality of the approximation depends mainly on the initialization. A random vector satisfying the constraints could be used. Instead, it seems to be more natural to compute a bipartite GED (bGED), i.e. to solve a LSAPE [START_REF] Bougleux | Linear sum assignment with edition[END_REF], and consider the corresponding optimal -assignment as initial vector in IPFP algorithm. In this way, results obtained by a bGED should be improved. This is analyzed experimentally in the following section.

IV. EXPERIMENTS

A. Chemical datasets

Table II shows results obtained by different state of the art methods to compute an estimation of the GED. These methods have been tested on four chemoinformatics datasets 1 , in the same way as in previous papers using these datasets [START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF], [START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF]. Each dataset is composed of molecular graphs. Costs have been set to 1 for any substitution operation, and 3 for insertions and deletions. It is important to notice that each graph adjacency matrix is randomly permuted in order to avoid some alignment bias which may provide approximations of the GED. The tested methods are:

• Exact GED based on A algorithm (line 1), restricted to very small graphs (not adapted to MAO and PAH). • Bipartite GED based on the Hungarian algorithm adapted to the LSAPE [START_REF] Bougleux | Linear sum assignment with edition[END_REF], with the substructure attached to each node defined by: incident edges and adjacent nodes [START_REF] Riesen | Bipartite graph matching for computing the edit distance of graphs[END_REF], [START_REF] Riesen | Approximate graph edit distance computation by means of bipartite graph matching[END_REF], a bag of random walks up to 3 edges (line 2) [START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF], a subgraph up to a radius 3 (line 3) [START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF]. • Approximate GED formalized as a QAP (line 4) [START_REF] Bougleux | A quadratic assignment formulation of the graph edit distance[END_REF] and as a QAPE (line 5), both solved by the IPFP algorithm with the same initialization provided by [START_REF] Gaüzère | Approximate Graph Edit Distance Guided by Bipartite Matching of Bags of Walks[END_REF], since initialization influences the approximation of the GED. Observe that these two last approaches obtain the best approximation results. Considering Alkane and Acyclic datasets, this is trivially inferred from average errors. For the two remaining datasets, we observe lower average edit distances which may correspond to better approximations since the exact GED is always over-estimated. In addition, we can note that computational times of quadratic approaches remain globally acceptable, especially compared to the bGED proposed in [START_REF] Carletti | Approximate graph edit distance computation combining bipartite matching and exact neighborhood substructure distance[END_REF] where most of the time is spent to compute the costs.

Moreover, quadratic approaches seem to have similar accuracies, while the approach presented in this paper improves clearly the computational time. This is analyzed through a second experiment consisting in 30 repetitions of each approach, for each dataset. As before, graph adjacency matrices are randomly permuted in order to avoid some alignment bias. Figure 1a shows the boxplot of the normalized average approximate GEDs. For each dataset, these average values are normalized by the maximal mean approximate GED given by max{max i { dQAPE i }, max i { dQAP i }} where di is the mean approximate GED obtained for the repetition i ∈ {1, . . . , 30}. Considering Acyclic and MAO datasets, the distribution are statistically similar according to a paired Student's test with a p-value equals to 0.01. Considering Alkane and PAH datasets, the difference is statistically significant but not so important for the approximation of GEDs (about 0.4%). From a computational time point of view, Figure 1b shows the normalized mean times according to a same initialization. We can clearly see that the proposed approach is faster than previous one, with a ratio around 2 for smaller graphs and 4 for larger graphs.

B. Synthetic datasets

To explore the behavior of the quadratic approaches when the number of nodes increases, we generate synthetic datasets having the same characteristics as MAO dataset but generalized to different graph sizes. These characteristics include node's and edge's labels distribution and ratio between number of edges and number of nodes. For a given number of nodes, each new synthetic graph is created according to these characteristics. This allows to retrieve characteristics of common molecular graphs. For each tested size, up to 100 nodes, a dataset composed of 100 synthetic graphs has been synthesized. Each of this graph, called source graph, is then associated to a target graph which has been obtained by removing one node and by substituting another one from the associated source graph. The exact GED between a source and a target is thus known by construction and is around 10.

Figure 2a shows the behavior of the quadratic approaches when the size increases. As we can see, the proposed approach based on a QAPE scales better in terms of computational times. Indeed, for graphs having 100 nodes, the computational ratio is over 5. Moreover, the average error is the same for both methods and remains quite stable on the set of sizes. It is worth noticing that the most important part of time is dedicated here to the initialization, realized as in the previous section. Figure 2b shows a similar result but with a fixed size for the first graph (20 nodes) and a varying size for the second one (from 20 to 500 nodes). The gain in computational times increases as the size of graphs also increases. Considering a second graph with 500 nodes, this gain corresponds to a computational time reduced by a factor 29. This is explained by the reduction of variable sizes from (n + m) 2 to (m + 1)(n + 1), and by the use of a LSAPE instead of a LSAP, at each iteration of IPFP algorithm (Sec. III-C).

V. CONCLUSION

Based on the notion of assignment with edition, this paper presents a binary quadratic programming problem whose global optimum is equal to the GED when graphs are simple, and when edit operations are restricted to substitutions, removals and insertions. The proposed quadratic program reformulates a previous expression of the GED as a QAP, by reducing the number of binary variables to (n + 1)(m + 1) and by relaxing the constraints associated to removals and insertions. This allows us to adapt an iterative algorithm, that initially finds an approximate solution to the QAP, to approximate the GED and reduce both space and time complexities. This is confirmed experimentally. Experiments also show that computational times are close to those of bGEDs, while providing distances closer to exact GEDs, on tested datasets.

Fig. 1 .

 1 Fig. 1. Comparison of IPFP algorithm for the QAP and for the QAPE.

Fig. 2 .

 2 Fig. 2. Total computational time in seconds taken by IPFP algorithm for the QAP (upper red line) and for the QAPE (bottom blue).

  Graph Edit Distance as a Quadratic Program

	Sébastien Bougleux	Benoit Gaüzère	Luc Brun
	Normandie Univ, UNICAEN	Normandie Univ, INSA Rouen	Normandie Univ, ENSICAEN
	CNRS, GREYC, France	LITIS, France	CNRS, GREYC, France

TABLE I EDIT

 I OPERATION INDUCED BY THE SIMULTANEOUS ASSIGNMENT OF

TABLE II ACCURACY

 II AND COMPUTATIONAL SCORES. d IS THE AVERAGE EDIT DISTANCE, e THE AVERAGE ERROR AND t THE AVERAGE COMPUTATIONAL TIME.

		Algorithm	d	Alkane e t	d	Acyclic e t	d	MAO t	d	PAH t
	1	A * 15		1.29	17		6.02
	2	[9], [10] 35	18	10 -3 35	18	10 -3 105	10 -3 138	10 -3
	3	[12] 33	18	10 -3 31	14	10 -2 49	10 -2 120	10 -2
	4	[13] 26	11	2.27	28	9	0.73	44	6.16	129	2.01
	5	QAP 20.2 4.7 0.002	20.8 3.4 0.004	33.9 0.019	53.9 0.038
	5 (this paper) QAPE 20.1 4.6	10 -3 20.8 3.5 0.002	33.9 0.005	53.6 0.010
	1
	0.996
	0.992
	0.988
	0.984
	(a) Repartition of normalized average mean edit distance for each repetition.

Datasets are available at https://iapr-tc15.greyc.fr/links.html


Nicolas Crouseilles
email: nicolas.crouseilles@inria.fr

Mohammed Lemou
email: mohammed.lemou@univ-rennes1.fr

Florian Méhats
email: florian.mehats@univ-rennes1.fr

Xiaofei Zhao

Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field

Keywords: Kinetic models, Vlasov-Poisson equation, Highly oscillatory, Four dimensional, Uniformly accurate, Particle-in-Cell, Two-scale formulation

In this work, we focus on the numerical resolution of the four dimensional phase space Vlasov-Poisson system subject to a strong external magnetic field. To do so, we consider a Particle-In-Cell based method, for which the characteristics are reformulated by means of the two-scale formalism, which is well-adapted to handle highly-oscillatory equations. Then, a numerical scheme is derived for the two-scale equations. The so-obtained scheme enjoys a uniform accuracy property, meaning that its accuracy does not depend on the small parameter. Several numerical results illustrate the capabilities of the method.

Introduction

In plasma devices like tokamaks, when the charged particles are subject to a strong external magnetic field, a drift phenomenon will happen in the plane orthogonal to the magnetic field direction. As a model to describe this phenomenon, we consider in this paper the long time Vlasov-Poisson equation with a strong external magnetic field in the four dimensional phase space. Namely, we consider a distribution function f ε (t, x, v) which depends on time t ≥ 0, space x = (x 1 , x 2 ) T ∈ R 2 and velocity v = (v 1 , v 2 ) T ∈ R 2 , which is solution to the following Vlasov-Poisson equation

∂ t f ε (t, x, v) + v ε • ∇ x f ε (t, x, v) + 1 ε E ε (t, x) + 1 ε v ⊥ • ∇ v f ε (t, x, v) = 0, (1.1a)
∇ x • E ε (t, x) = ρ ε (t, x) -n i , ρ ε (t, x) := R 2 f ε (t, x, v)dv, (1.1b)
f ε (0, x, v) = f 0 (x, v), x, v ∈ R 2 , (1.1c)
where E ε is the electric field solution of the Poisson equation (1.1b), and where f 0 is a given initial data.

We denote by v ⊥ as

v ⊥ = Jv = (v 2 , -v 1 ) T , with J = 0 1 -1 0 .
The real-valued scalar function f ε (t, x, v) is the distribution function of electrons. The given parameters ε > 0 is inversely proportional to the strength of the magnetic field and n i > 0 denotes the given ion density of the background which ensures global quasi-neutrality and hence the solvability of the Poisson equation (1.1b).

When the plasma is subject to a strong magnetic field, i.e. the parameter ε becomes very small, the solution f ε (t, x, v) of the Vlasov-Poisson equation (1.1) exhibits fast oscillations in time with wavelength O(ε 2 ). On the other side, there exists a slower dynamics of the so-called guiding-center, which is essentially given by the "E-cross-B" drift velocity (see [START_REF] Lee | Gyrokinetic approach in particle simulation[END_REF][START_REF] Dubin | Nonlinear gyrokinetic equations[END_REF]). Thus, the system (1.1) is strongly multi-scale which implies strong restrictions on the numerical parameters when one wants to simulate the solution of (1.1).

When one is interested in the regime ε < < 1, the asymptotic model is sufficient to describe the system in this case; and since it does not contain any stiff term, standard numerical methods can be employed to approximate it. The solution f ε of (1.1) has been proved to converge in a suitable sense towards f when ε goes to zero (see [START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF]) which solves

∂ t f (t, x, v) + E ⊥ (t, x) • ∇ x f (t, x, v) + 1 2 ∇ x • E(t, x)v ⊥ • ∇ v f (t, x, v) = 0, (1.2a)
∇ x • E(t, x) = R 2 f (t, x, v)dv -n i . (1.2b) (1.2c)
As usual, the electric field E derives from a potential φ: E = -∇ x φ. By integrating (1.2a) with respect to v ∈ R 2 , we further get the well-known guiding-center model

∂ t ρ(t, x) + E ⊥ (t, x) • ∇ x ρ(t, x) = 0, t > 0, x ∈ R 2 , (1.3a) ∇ x • E(t, x) = ρ(t, x) -n i , (1.3b)
ρ(0, x) = R 2 f 0 (x, v)dv. (1.3c)
The unknown ρ is the limit of the density ρ ε as ε → 0 (see [START_REF] Saint-Raymond | The gyro-kinetic approximation for the Vlasov-Poisson system[END_REF][START_REF] Bostan | The Vlasov-Maxwell system with strong initial magnetic field. Guiding-center approximation[END_REF] for mathematical justification of this limit). As said before, these models are only valid in the asymptotic regime ε < < 1.

From a numerical point of view, the passage from (1.1) to (1.3) has been explored by recent works, within the framework of Particle-In-Cell (PIC) method. The PIC framework (see [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Jacobs | Implicit-Explicit time integration of a high-order particle-in-cell method with hyperbolic divergence cleaning[END_REF]) enables to focus on the construction of a numerical integrator in time of the characteristics equations of (1.1). In particular, in [START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF], an exponential time integrator has been proposed whereas in [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF], the authors propose a high order IMEX type scheme. Both approaches are uniformly stable with respect to ε, and degenerate when ε goes to zero to a consistent scheme with the asymptotic model (1.3). In particular, the numerical parameters can be chosen independently from ε. However, this so-called Asymptotic Preserving property (see [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF]) does not guarantee the good behavior of the numerical scheme in the intermediate regime. The main goal of the present work is to construct a Uniformly Accurate (UA) numerical scheme with respect to ε. More precisely, in this work we will also use a Particle-In-Cell method, in which the characteristics equation will be reformulated using the two-scale formalism (as introduced in [START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF]). In addition to the uniform stability, this strategy enables to capture the highly oscillatory (in time) behavior of the solution, and to derive a uniformly accurate Particle-In-Cell method for ε ∈ (0, 1], which means that the accuracy of the scheme does not depend on ε.

The strategy we will use in this work relies on augmented characteristics equations, where the two time scales t/ε 2 and t are separated. The solution of the original characteristics equation are recovered from the two-scale solution by evaluating it on the diagonal (t, τ = t/ε 2 ). Although the overall strategy follows the one introduced in [START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF], several improvements have to be done to handle the problem (1.1). First, the model (1.1) under study involves two stiff terms: some transport terms (in x, v) of magnitude 1/ε and a transport term in v of magnitude 1/ε 2 . This contrasts with the previous study where only one stiffness of magnitude 1/ε was involved, and naturally induces several difficulties. For instance, as we shall see, even after the filtering step, a nonlinear stiff term is still present in the reformulated equations. As this stiff term strongly depends on the self-consistent electric field, which is computed numerically, the numerical error is amplified by the stiffness factor 1/ε and the uniform accuracy is difficult to obtain (contrary to what has been done in [START_REF] Lemou | Uniformly accurate numerical schemes for the nonlinear Dirac equation in the nonrelativistic limit regime[END_REF] where the analogous term was linear and then was considered implicit in time). To overcome this difficulty, a suitable change of unknown is performed, which enables to reformulate the two-scale model in order to fit with the theoretical framework of [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF]. A final and important step is the suitable preparation of the initial data for the two-scale model. As proved in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF], a suitable choice of the initial condition (in the additional variable τ ) ensures that the time derivatives of the two-scale model are uniformly bounded with respect to ε. In the present work, a second order exponential integrator is used for the two-scale model, for which the local truncation error only depends on a second time derivative of the solution (whereas the second order scheme used in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF] induces the third time derivative). Then, a first order preparation of the initial data is sufficient to ensure its uniform boundedness and as a consequence, the uniform accuracy of the numerical scheme is obtained. Note that this first order correction should take into account the dependence in τ of the electric field and a tricky expansion with respect to ε of the self-consistent electric field is then required.

It is worth mentionning that this UA PIC scheme is explicit in time and does not require any iterative algorithm. Moreover, as a byproduct of the uniform accuracy, the scheme has the good asymptotic behavior when ε << 1, since it is consistent with the asymptotic model (1.2). Several numerical tests are performed to validate and illustrate the capabilities of the UA PIC scheme in the 4 dimensional phase space setting. In particular, a few number of points in the additional τ direction are needed to produce accurate numerical results. Note that this τ direction is easily parallelized using MPI paradigm and the extra cost turns out to be negligible compared to the cost of a standard Particle-In-Cell method.

The rest of the paper is organized as follows. In the next Section, the UA PIC method is presented in details after a brief introduction of PIC methods. Then, some numerical results are given in Section 3 and conclusions are drawn in Section 4.

Uniformly accurate particle-in-cell method

In this section, we shall begin by briefly presenting the framework of the Particle-in-Cell (PIC) method where more details could be found in [START_REF] Sonnendrücker | Numerical Methods for Vlasov-Maxwell Equations[END_REF][START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF], and then we shall introduce our uniformly accurate (UA) integrator towards the characteristics equation.

Framework of PIC

The starting point of the Particle-in-Cell method is to approximate the unknown distribution f ε (t, x, v) of (1.1) by a sum of Dirac masses centred at (x k (t), v k (t)) with weights ω k > 0 for k = 1, . . . , N p and N p ∈ N as

f ε p (t, x, v) = Np k=1 ω k δ(x -x k (t))δ(v -v k (t)), t ≥ 0, x, v ∈ R 2 . (2.1)
Inserting (2.1) into (1.1), in the sense of distribution, for k = 1, . . . , N p , each particle obeys the characteristics equation

ẋk (t) = v k (t) ε , (2.2a
)

vk (t) = E ε (t, x k (t)) ε + v ⊥ k (t) ε 2 , t > 0, (2.2b) x k (0) = x k,0 , v k (0) = v k,0 . (2.2c)
The weight ω k and initial values of the particles x k,0 , v k,0 for k = 1, . . . , N p are prescribed according to the initial distribution f 0 (x, v) given by (1.1c). To determine the weight, by integrating (2.1) at t = 0 in whole space we require

Np k=1 ω k = R 2 ×R 2 f 0 (x, v)dxdv.
Thus, a simple choice of uniform weight for all particles would be

ω k = 1 N p R 2 ×R 2 f 0 (x, v)dxdv, k = 1, . . . , N p .
The initialization of x k,0 , v k,0 for k = 1, . . . , N p is a rather well-established and classical sampling issue. It could be done by the Monte Carlo type rejection sampling method as a general approach. The detailed process can be found in standard statistics textbooks. As a consequence, this approach introduces a noise of order 1/ N p (see [START_REF] Sonnendrücker | Numerical Methods for Vlasov-Maxwell Equations[END_REF]). Note however that when f 0 is of variable separation form, i.e.

f 0 (x, v) = χ 1 (x)χ 2 (v)
where χ 2 is a Gaussian, specific pseudo-random methods can be used. We omit the details here for brevity and refer the reader to [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF]. The characteristics equations for particles (2.2) are coupled to the Poisson equation (1.1b) through the electrical field E ε . Once the positions {x k (t)} k=1,...,Np of the particles are obtained at time t > 0, one needs to evaluate the approximated density

ρ ε p (t, x) = Np k=1 ω k δ(x -x k (t)) ≈ ρ ε (t, x) x ∈ R 2 , (2.3)
then solve the Poisson equation

∇ x • E ε p (t, x) = ρ ε p (t, x) -n i , x ∈ R 2
, on a mesh grid of the spatial domain, with step size ∆x. This provides an approximation E ε p (t, x) ≈ E ε (t, x). Finally, the values of E ε p (t, x k (t)) are interpolated at each particle position. In practice, the Dirac delta function δ(x) is replaced by smooth functions S(x) (particle shape function). For example in 1D, it is approximated by the B-spline function S m (x) of order m ∈ N [START_REF] Sonnendrücker | Numerical Methods for Vlasov-Maxwell Equations[END_REF]:

S 0 (x) :=    1 ∆x , |x| ≤ ∆x 2 , 0, else, , S m (x)
:= 1 ∆x

x+ ∆x 2 x-∆x 2 S m-1 (y)dy, m ≥ 1. (2.4)
The case in two dimensions is done by tensor product. The B-spline function S m is defined locally in space which is preferred from numerical point of view, but globally it is only a C m-1 function for m ≥ 1. A smooth but global basis has been considered in [START_REF] Jacobs | Implicit-Explicit time integration of a high-order particle-in-cell method with hyperbolic divergence cleaning[END_REF] to get high order accuracy. The classical PIC scheme is hence completed by a particle pusher for (2.2).

As is known, either the classical numerical methods or the asymptotic preserving methods would be problematic for solving oscillatory systems as (2.2) when ε is in a certain regime, which in turns make the PIC scheme not efficient or less accurate for ε ∈ (0, 1]. Thus, the main difficulty here is to solve the characteristic equation by some efficient integrator that works uniformly well for all 0 < ε ≤ 1 and embed the integrator harmoniously into the PIC framework. We are going to apply the two-scale formulation approach developed in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF] to derive a uniformly accurate PIC method.

Two-scale formulation for characteristics

This part is devoted to the derivation of the two-scale formulation of the highly-oscillatory ODE system (2.2). To do so, we follow two main steps: after a classical filtering step, the two-scale system can be obtained and then, a suitable change of unknown enables to recast the so-obtained two-scale system in a formulation like in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF]. Note that for simplicity of notation, in this subsection we shall omit the subscript k for particles in (2.2).

First, we introduce a new variable y in order to filter out the main oscillation as

y(t) = e -tJ/ε 2 v(t), t ≥ 0, with e -tJ = cos t -sin t sin t cos t , (2.5)
so that the characteristic equation (2.2) becomes

ẋ(t) = 1 ε e tJ/ε 2 y(t), (2.6a) ẏ(t) = 1 ε e -tJ/ε 2 E(t, t/ε 2 , x(t)), (2.6b
)

x(0) = x 0 , y(0) = v 0 . (2.6c)
Then, we derive a suitable two-scale formulation by separating the fast scale τ = t/ε 2 from the slow one t, as already done in [START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Lemou | Uniformly accurate numerical schemes for the nonlinear Dirac equation in the nonrelativistic limit regime[END_REF]. We consider an augmented system of equations satisfied by the augmented solution (X k (t, τ ), Y k (t, τ )) (k = 1, . . . , N p ) which is periodic with respect to τ (of period 2π), and which coincides with the original solution (x k (t),

y k (t)) (k = 1, . . . , N p ) on the diagonal τ = t/ε 2 x k (t) = X k t, t ε 2 , y k (t) = Y k t, t ε 2 , t ≥ 0, k = 1, . . . , N p .
The two-scale equation satisfed by

(X k (t, τ ), Y k (t, τ )) reads (removing the subscript k) ∂ t X(t, τ ) + 1 ε 2 ∂ τ X(t, τ ) = 1 ε e τ J Y (t, τ ), (2.7a) ∂ t Y (t, τ ) + 1 ε 2 ∂ τ Y (t, τ ) = 1 ε e -τ J E(t, τ, X(t, τ )), (2.7b)
where the electric field is given by

∇ x • E(t, τ, x) ≈ Np =1 ω δ(x -X (t, τ )) -n i . (2.8)
The augmented equation (2.7) are to be supplemented with the following initial data

X(0, 0) = x(0), Y (0, 0) = y(0).
Note that we can recover the original solution of (2.2) as

x(t) = X t, t ε 2 , v(t) = e tJ/ε 2 Y t, t ε 2 , t ≥ 0.
Solving (2.7) requires to solve numerically the augmented Poisson equation (2.8). Then, at each time, an approximation error of E(t, τ, x) will be introduced in (2.7). This error (which depends on the number of particles N p and the mesh size ∆x in x-space) is amplified by the factor 1/ε on the right hand side of (2.7b), due to the diffusion scaling. As a matter of fact, the extra error apart from the time discretization will make the scheme non-uniformly accurate for ε ∈ (0, 1]. Note that this lack of accuracy does not occur when E is given or explicitely depends on the solution. This later case corresponds to [START_REF] Lemou | Uniformly accurate numerical schemes for the nonlinear Dirac equation in the nonrelativistic limit regime[END_REF].

To overcome this essential difficulty, we seek for a suitable reformulation of (2.7). Let us introduce two new unknown functions U ± (t, τ ):

U + (t, τ ) := X(t, τ ) + εJe τ J Y (t, τ ), U -(t, τ ) := -εJY (t, τ ), t ≥ 0, τ ∈ [0, 2π]. (2.9)
The unknown (X(t, τ ), Y (t, τ )) is obtained by the inverse change of variables

X(t, τ ) = U + (t, τ ) + e τ J U -(t, τ ) , Y (t, τ ) = 1 ε JU -(t, τ ). (2.10)
Combining (2.9) with (2.7) leads to

∂ t U ± + 1 ε 2 ∂ τ U ± = F ± (t, τ, U + , U -), (2.11)
with

F + (t, τ, U + , U -) := JE t, τ, U + + e τ J U -, F -(t, τ, U + , U -) := -Je -τ J E t, τ, U + + e τ J U -, (2.12)
with the initial condition

U + (0, 0) = x 0 + εJv 0 =: u 0 + , U -(0, 0) = -εJv 0 =: u 0 -, (2.13)
x 0 , v 0 being the initial condition of the original problem (2.2). Note that E is computed from (U ± ) , = 1, . . . , N p by the Poisson equation

∇ x • E(t, τ, x) = Np =1 ω δ x -(U + + e τ J U -) . (2.14)
Again, from the solution of (2.11), the solution of (2.2) is recovered via the following relations

x(t) = U + (t, t/ε 2 ) + e tJ/ε 2 U -(t, t/ε 2 ), v(t) = 1 ε e tJ/ε 2 JU -(t, t/ε 2 ). (2.15)
Comparing to (2.7), now (2.11) has bounded coefficients in front of the nonlinearity and it enters in the framework studied [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF]. In [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF], it is proven that, the initial condition U ± (t = 0, τ ) has to be chosen under a suitable expansion to ensure that the time derivatives of U ± are uniformly bounded with respect to ε. This is of importance when one applies a numerical discretization to (2.11) since a uniform accuracy property can then proved for the numerical scheme (see [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Lemou | Uniformly accurate numerical schemes for the nonlinear Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Nonlinear Geometric Optics method based multi-scale numerical schemes for highlyoscillatory transport equations[END_REF]). The construction of the suitable initial data U ± (t = 0, τ ) is based on a Chapman-Enskog type expansion and is detailed in the next subsection.

Derivation of the initial condition

Following the strategy proposed in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Lemou | Uniformly accurate numerical schemes for the nonlinear Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Nonlinear Geometric Optics method based multi-scale numerical schemes for highlyoscillatory transport equations[END_REF], we will derive a suitable initial condition for (2.11). To do so, we will perform a Chapman-Enskog expansion of U ± (t, τ ) so that U ± (0, 0) = u 0 ± .

Chapman-Enskog expansion of U ± . We denote Lu(τ ) = ∂ τ u(τ ) for some periodic function u(τ ) on [0, 2π] and introduce the average operator Π as

Πu = 1 2π 2π 0 u(s)ds.
Then L is invertible on the set of periodic functions with zero average, i.e. for u(τ ) with Πu = 0,

L -1 u(τ ) = (I -Π) τ 0 u(s)ds.
We introduce the decomposition

U ± (t, τ ) = U ± (t) + h ± (t, τ ), with U ± (t) = ΠU ± (t, τ ),
then we have the equation for the macro part of the solution

∂ t U ± (t) = ΠF ± (t, τ, U + (t) + h + (t, τ ), U -(t) + h -(t, τ )),
and the equation for the micro part

∂ t h ± (t, τ ) + 1 ε 2 Lh ± (t, τ ) = (I -Π)F ± (t, τ, U + (t) + h + (t, τ ), U -(t) + h -(t, τ )). (2.16)
By inverting L in the micro part equation (2.16), we get

h ± (t, τ ) = ε 2 AF ± (t, τ, U + (t) + h + (t, τ ), U -(t) + h -(t, τ )) -ε 2 L -1 ∂ t h ± (t, τ ),
where

A := L -1 (I -Π). By assuming ∂ t h ± = O(1), we get h ± = O(ε 2 ). By further assuming ∂ 2 t h ± = O(1)
, we get a fourth order expansion of U ± as

U ± (t, τ ) = U ± (t) + ε 2 AF ± (t, τ, U + (t), U -(t)) + O(ε 4 ), (2.17)
where we recall that F ± (t, τ, U + (t), U -(t)) are given by

F + (t, τ, U + (t), U -(t)) = JE t, τ, U + + e τ J U + , (2.18)
F -(t, τ, U + (t), U -(t)) = -Je -τ J E t, τ, U + + e τ J U + . (2.19)
One has to point out that the electric field E is computed through the Poisson equation (2.14), and then evaluated at (U + + e τ J U + ). To get an expansion of U ± as a function of U ± , one then needs to also expand the electric field E solution of (2.14). To do so, we approximate U ± in the right hand side of (2.14) by U ± . We will see that this approximation provides a second order approximation of E, as required from the expansion (2.17).

Denoting by X 1st the following quantity

X 1st (t, τ ) := U + (t) + e τ J U -(t), (2.20)
one has X = X 1st + O(ε 2 ). Indeed, on the one side, a Chapman-Enskog expansion of X solution of (2.7) reads

X(t, τ ) = X(t) + εL -1 e τ J Y (t, τ ) + O(ε 2 ) = X(t) -εJe τ J Y (t) + O(ε 2 ).
On the other side, one can express X 1st defined in (2.20) in terms of X and Y using (2.9)

X 1st (t, τ ) = X(t) + εJΠ(e τ J Y (t, τ )) -εe τ J JY (t) = X(t) -εJe τ J Y (t) + O(ε 2 ).
We deduce from the two last equation that X = X 1st + O(ε 2 ). Then, using the Poisson equation (2.8), one can compute

E 1st ∇ x • E 1st (t, τ, x) ≈ Np k=1 ω k δ(x -X 1st k (t, τ )) -n i , (2.21)
and we deduce E(t, τ, x) = E 1st (t, τ, x) + O(ε 2 ). Finally E 1st enables to get F 1st ± as in (2.18)-(2.19):

F 1st + (t, τ, U + (t), U -(t)) = JE 1st t, τ, U + (t) + e τ J U -(t) , (2.22)
F 1st -(t, τ, U + (t), U -(t)) = -Je -τ J E 1st t, τ, U + (t) + e τ J U -(t) . (2.23)
We end up by defining U 2nd as

U 2nd ± (t, τ ) = U ± (t) + ε 2 AF 1st ± (t, τ, U + (t), U -(t)). ( 2

.24)

Link with the original initial condition u 0 ± . Now, we need to compute U ± so that U 2nd ± (0, 0) = u 0 ± where U 2nd is defined by (2.24). Note that u 0 ± is defined by (2.13) as a function of (x 0 , v 0 ) which is the initial data of the original problem. Evaluating (2.24) at t = τ = 0 and using U 2nd ± (0, 0) = u 0 ± , one has directly

u 0 ± -U ± (t = 0) = O(ε 2
). Then, we can rewrite (2.24) at t = 0 as

U 2nd ± (0, τ ) = U ± (0) + h 1st ± (τ ) + O(ε 4 ), with h 1st ± (τ ) = ε 2 AF 1st ± (0, τ, u 0 + , u 0 -), (2.25)
and F 1st ± defined by (2.22)-(2.23). Using again the condition U 2nd ± (0, 0) = u 0 ± , we get from (2.25) at τ = 0:

U ± (0) = u 0 ± -h 1st ± (0)
. Replacing U ± (0) by this latter expression in (2.25) leads to

U 2nd ± (0, τ ) = u 0 ± -h 1st ± (0) + h 1st ± (τ ), with h 1st ± (τ ) := ε 2 AF 1st ± (0, τ, u 0 + , u 0 -), (2.26) with F 1st + (0, τ, u 0 + , u 0 -) = JE 1st 0, τ, u 0 + + e τ J u 0 -, F 1st -(0, τ, u 0 + , u 0 -) = -Je -τ J E 1st 0, τ, u 0 + + e τ J u 0 -,
where E 1st is computed with the following Poisson equation

∇ x • E 1st (0, τ, x) = Np k=1 ω k δ x -(u 0 + ) k + e τ J (u 0 -) k -n i ,
and where u 0 + = x 0 + εJv 0 , u 0 -= -εJv 0 . We conclude this part with the following proposition which ensures that the solution of the two-scale equations (2.11) with the suitable initial condition (2.26) is smooth with respect to ε. Proposition 2.1. Let us consider the following problem satisfied by

U ± (t, τ ), t ∈ [0, T ], T > 0, τ ∈ [0, 2π] ∂ t U ± + 1 ε 2 ∂ τ U ± = F ± (t, τ, U + , U -),
with the initial condition defined by (2.26) and the above notations. Then, there exists a constant C > 0 independent of ε such that,

∂ k t U ± L ∞ t,τ ≤ C, k = 0, 1, 2.
(2.27)

The proof will not be done here but can be obtained in a analogous way as in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF].

An exponential integrator

Now we present a numerical scheme for the discretization of the characteristics (2.11). We shall integrate the two-scale equation (2.11) by an exponential integrator.

Let ∆t > 0 be the time step and denote t n = n∆t for n ≥ 0. Discretizating the τ -direction as τ j = j∆τ, j = 0, 1, . . . , N τ , with ∆τ = 2π/N τ and N τ some positive even integer, and applying the Fourier transform to (2.11) on [0, 2π], we get for = -N τ /2, . . . , N τ /2 -1,

d dt (U ± ) (t) + i ε 2 (U ± ) (t) = (F ± ) (t), t > 0, (2.28)
where

U ± (t, τ ) = Nτ /2-1 =-Nτ /2 (U ± ) (t)e i τ , F ± (t, τ, U + , U -) = Nτ /2-1 =-Nτ /2 (F ± ) (t)e i τ .
Integrating (2.28) from t n to t n+1 (n ≥ 0), we get

(U ± ) (t n+1 ) = e -i ∆t ε 2 (U ± ) (t n ) + ∆t 0 e -i ε 2 (∆t-s) (F ± ) (t n + s)ds. (2.29)
We approximate the above integral by the following quadrature for n ≥ 1,

(U ± ) (t n+1 ) ≈ e -i ∆t ε 2 (U ± ) (t n ) + ∆t 0 e -i ε 2 (∆t-s) (F ± ) (t n ) + s d dt (F ± ) (t n ) ds ≈ e -i ∆t ε 2 (U ± ) (t n ) + p (F ± ) (t n ) + q 1 ∆t (F ± ) (t n ) -(F ± ) (t n-1 ) ,
where

p := ∆t 0 e -i ε 2 (∆t-s) ds =    iε 2 e -i ∆t ε 2 -1 , = 0, ∆t, = 0, q := ∆t 0 e -i ε 2 (∆t-s) sds =        ε 2 2 ε 2 -ε 2 e -i ∆t ε 2 -i ∆t , = 0, ∆t 2 2 , = 0.
For the first iteration, the starting value (U ± ) (t 1 ) (with t 1 = ∆t) is computed using a prediction-correction scheme. The first order prediction writes

(U ± ) * (t 1 ) := e -i ∆t ε 2 (U ± ) (0) + p (F ± ) (0),
which gives a prediction of the nonlinearity at t 1 = ∆t as

F ± (t 1 , τ, U * + (t 1 , τ ), U * -(t 1 , τ )) = Nτ /2-1 =-Nτ /2 (F ± ) * (t 1 )e i τ .
Then, we do a correction step

(U ± ) (t 1 ) ≈ e -i ∆t ε 2 (U ± ) (0) + p (F ± ) (0) + q 1 ∆t (F ± ) * (t 1 ) -(F ± ) (0) .
At each time level t n (n ≥ 0), when the approximated U ± (t n , τ ) is obtained as above, we get

X(t n , τ ) = U + (t n , τ ) + e τ J U -(t n , τ ), t ≥ 0, τ ∈ [0, 2π].
and we compute

E(t n , τ, x) = -∇ x Φ ε (t n , τ, x), -∆Φ ε (t n , τ, x) ≈ Np k=1 ω k δ(x -X k (t n , τ )) -n i ,
then the evaluation of nonlinearity F ± follows from (2.12)

F + (t n , τ, U + (t n , τ ), U -(t n , τ )) = JE t n , τ, U + (t n , τ ) + e τ J U -(t n , τ ) , F -(t n , τ, U + (t n , τ ), U -(t n , τ )) = -Je -τ J E t n , τ, U + (t n , τ ) + e τ J U -(t n , τ ) .
In all, the detailed numerical scheme of the solver towards (2.11) reads as follows. Denote (U ± ) n (τ ) ≈ U ± (t n , τ ) for n ≥ 0 and let U 0 ± = U ± (0, τ ) the initial condition given by (2.26). We update the U n ± for n ≥ 1 as

(U ± ) 1 = e -i ∆t ε 2 (U ± ) 0 + p (F ± ) 0 + q 1 ∆t (F ± ) * ,1 -(F ± ) 0 , (2.30a)
(U ± ) n+1 = e -i ∆t ε 2 (U ± ) n + p (F ± ) n + q 1 ∆t (F ± ) n -(F ± ) n-1 , ∀n ≥ 1, (2.30b)
where

U n ± (τ ) = Nτ /2-1 =-Nτ /2 (U ± ) n e i τ , F n ± (τ ) = Nτ /2-1 =-Nτ /2 (F ± ) n e i τ , F * ,1 ± (τ ) = Nτ /2-1 =-Nτ /2 (F ± ) * ,1
e i τ , n ≥ 0, and

F n + (τ ) = JE n τ, U n + (τ ) + e τ J U n -(τ ) , F n -(τ ) := -Je -τ J E n τ, U n + (τ ) + e τ J U n -(τ ) , n ≥ 0, F * ,1 + (τ ) = JE * ,1 τ, U * ,1 + (τ ) + e τ J U * ,1 -(τ ) , F * ,1 -(τ ) := -Je -τ J E 1, * τ, U * ,1 + (τ ) + e τ J U * ,1 -(τ ) , with (U ± ) * ,1 = e -i ∆t ε 2 (U ± ) 0 + p (F ± ) 0 , U * ,1 ± (τ ) = Nτ /2-1 =-Nτ /2 (U ± ) * ,1 e i τ .
The numerical electrical field E n is given by collecting all the two-scale characteristics (U n ± ) k , k = 1, . . . , N p and solving the Poisson equation:

∇ x • E n (τ, x) = Np k=1 ω k S(x -X n k (τ )) -n i , X n k (τ ) = U n +,k (τ ) + e τ J U n -,k (τ ), ∇ x • E * ,1 (τ, x) = Np k=1 ω k S(x -X * k (τ )) -n i , X * k (τ ) = U * ,1 +,k (τ ) + e τ J U * ,1 -,k (τ ), (2.31)
where S(x) is the regularised basis (the choice of quintic splines has been done in our computations).

We conclude this section with the following proposition which states that the so-obtained numerical scheme is uniformly accurate with respect to ε. Proposition 2.2. Let U ± the solution of (2.11) on [0, T ] × [0, 2π] subject to the initial condition (2.26). Let (U n ± ), n ≥ 0 defined by the numerical scheme (2.30) with the initial data (2.26). Then, the following estimate holds

sup ε∈]0,1] U ± (t n , •) -U n ± (•) L ∞ τ ≤ C∆t 2 ,
with C > 0, for all n = 0, . . . , N with N ∆t = T is the final time.

Proof. By Taylor expansion (with integral remainder) of V (t n + s) := (F ± ) (t n + s) in (2.29), we get

V (t n + s) = V (t n ) + s ∆t [V (t n ) -V (t n-1 )] + 1 0 (1 -t)s 2 V (t n + st)dt + s 1 0 (1 -t)∆tV (t n -t ∆t)dt.
Then, the exact solution (U ± )(t) satisfies the scheme (2.30) up to the following remainder term R R := ∆t 0 e -i ε 2 (∆t-s) s 2 1 0

(1 -t)V (t n + st)dtds

+ ∆t 0 e -i ε 2 (∆t-s) s 1 0 ∆t(1 -t)V (t n -t ∆t)dtds.
Thanks to Proposition 2.1, ∂ 2 t U L ∞ t,τ is uniformly bounded (the derivatives of F are bounded, as in [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF]) so that

|R| ≤ ∆t 0 s 2 ds + ∆t ∆t 0 sds V L ∞ t,τ ≤ C∆t 3 .

Now, we define the error E

n := | (U ± )(t n ) -(U ± ) n
|, which satisfies the following estimate

E n+1 ≤ E n + C ∆t 0 (E n + E n-1 )ds + C∆t 3 ≤ E n (1 + C∆t) + C∆tE n-1 + C∆t 3 . (2.32)
We proceed by induction to prove there exisits C > 0 such that E n ≤ C∆t 2 (1 + C∆t) n , ∀n ≥ 1. For the first iteration, from E 0 = 0 then we get E 1 ≤ C∆t 2 (since a prediction-correction scheme is used for the first iteration). Thanks to this estimate, we conclude the proof using the following

E n ≤ C∆t 2 (1 + C∆t) n ≤ C∆t 2 e Cn∆t ≤ C∆t 2 e CT ,
where T = N ∆t is the final time.

The proposed exponential integrator (2.30)-(2.31) is self-consistent for solving (2.11). During the computation, the scheme runs with U n ± (τ ) only. It does not require approximations of the original components x(t), y(t) for (2.6) or X(t, τ ), Y (t, τ ) for (2.7). In addition to the second order uniform accuracy in time and the spectral accuracy in τ , the scheme (2.30)-(2.31) is fully explicit and efficient thanks to fast Fourier transform in τ and the locality of the regularised basis S(x).

With (U n ± ) k , k = 1, . . . , N p from the exponential integrator (2.30)-( 2.31), one can obtain the approximation of the solution of the original characteristics (2.2) as

x k (t n ) ≈ U n + (t n /ε 2 ) + e Jtn/ε 2 U n -(t n /ε 2 ) k = X n k t n ε 2 , v k (t n ) ≈ J ε e J tn ε 2 U n -,k t n ε 2 , n ≥ 1.
Then together with the PIC approximation (2.1), we complete the numerical scheme of a two-scale PIC method for solving the Vlasov-Poisson equation (1.1).

Numerical results

In this section, we first test the accuracy of the proposed UA PIC method (2.30)-(2.31) and then conduct some numerical experiments to explore the dynamics of the Vlasov-Poisson equation (1.1) in the limit regime.

Example 1

We consider the Vlasov-Poisson equation (1.1) with the Kelvin-Helmholtz instability type initial data [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Shoucri | A two-level implicit scheme for the numerical solution of the linearized vorticity equation[END_REF][START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF],

f 0 (x, v) = 1 2π (1 + sin(x 2 ) + η cos(kx 1 )) e -|v| 2 2 , (3.1)
on a computational domain for x = (x 1 , x 2 ) as Ω = [0, 2π/k] × [0, 2π] for some k, η > 0. Periodic boundary conditions are assumed in x.

Accuracy test. In order to show the importance of the preparation of the initial data for the UA property, we consider the following initial data for each two-scale characteristic (2.11):

• U ± (0, τ ) = u 0 ± (we will refer as U 0th );

• U ± (0, τ ) = U 1st ± ; here the electric field E 1st is replaced in the Chapman-Enskog expansion (2.26) by the initial electric field

E(t = 0, x) solution of ∇ x • E(t = 0, x) = Np =1 ω δ(x -x ,0 ); • U ± (0, τ ) = U 2nd ±
given by (2.26).

We shall test the performance of the proposed scheme (2.30)-(2.31) under the two inefficient choices of initial data U 0th ± and U 1st ± . We shall compare the results with the one obtained from the two-scale PIC (2.30)-(2.31) under the well-prepared initial data U 2nd ± (2.26). We choose η = 0.05, k = 0.5 in (2.26) for this accuracy test. We discretize the spatial domain Ω with 64 points in x 1 -direction and 32 points in x 2 -direction. The Poisson equation (1.1b) is solved on Ω with periodic boundary conditions by the fast Fourier transform where the background density n i is such that the average in space of total charge vanishes. We take the basis function S(x) for the PIC method as the fifth order B-spline (same in other numerical experiments), i.e. m = 5 in (2.4). The high order spline function is chosen here to reach a high accuracy in spatial approximation so that the temporal error can be distinguished. We choose the numerical grid points in τ -direction as N τ = 32. We take N p = 204800 and generate the initial position and velocity (x k,0 , v k,0 ) of each particle by means of the rejection sampling.

We compute the error on the density function ρ ε (t, x) = R 2 f ε (t, x, v)dv, x ∈ Ω at the final time t = 1 where a reference solution is obtained numerically via the UA PIC with very small time step, e.g. ∆t = 10 -3 . We measure the error under maximum norm in space. The temporal error test results of the proposed twoscale PIC method (2.30)-(2.31) with the three types of initial data U 0th ± , U 1st ± , U 2nd ± are given in Figure 1. We can observe the influence of the initial data on the convergence of the method. When the initial data is chosen as U 2nd ± , the uniform accuracy is well recovered whereas order reduction (or even no convergence) can be observed when the initial data is not well-prepared, in particular in the intermediate regime. We also study the numerical error to further address the error from the velocity part, considering the following quantities

ρ ε v,1 (t, x) := R 2 (|v 1 | + |v 2 |)f ε (t, x, v)dv, ρ ε v,2 (t, x) := R 2 |v| 2 f ε (t, x, v)dv, x ∈ Ω.
The errors of the UA PIC regarding ρ ε v,1 and ρ ε v,2 at t = 1 are shown in Figure 2, with the initial condition U 2nd ± (under the similar setup as previously). The uniform accuracy is still observed on these quantities. Then, in the following, UA PIC scheme will always refer to the scheme (2.30)-(2.31) with the initial data U 2nd ± . Figure 3 investigates the error in τ -direction (using the maximum norm) by comparing the solution obtained by the UA PIC scheme under ∆t = 10 -3 and different N τ to a reference one obtained under N τ = 32 (the other numerical parameters are the same as previously). Figure 3 illustrates the spectral accuracy of the scheme in the τ -direction, meaning that a few number of grid points are sufficient to obtain Finally, the time history of the relative total energy (H(t) -H(0))/H(0) is displayed in Figure 4, where H(t) is a conserved quantity defined by

∆ t=0.2 ∆ t=0.2/2 ∆ t=0.2/2 2 ∆ t=0.2/2 3 ∆ t=0.2/2 4
H(t) = 1 2 R 2 Ω |v| 2 f ε (t, x, v)dxdv + 1 2 Ω |E ε (t, x)| 2 dx.
The error of the total energy, i.e. |H(t n ) -H(0)|/H(0) during the computation of the UA PIC (i.e. with U 2nd ± ) is given in Figure 4 under time step ∆t = 0.1 and ∆t = 0.05 for several ε. We can observe that the error is well preserved and in particular, its amplitude is divided by 4 when the time step is divided by 2, as expected.

Numerical exploration. Next, we carry out some numerical explorations on the Vlasov-Poisson equation (1.1). We firstly investigate numerically the convergence of the Vlasov-Poisson equation (1.1) to its limit model (1.2). We work under the initial data (3.1) with η = 0.05, k = 0.5. We solve the asymptotic model (1.2) satisfied by f (t, x, v) using a standard PIC method (a leap-frog scheme is used for the characteristics).

We study the following quantities

ρ ε (t, x) -ρ(t, x) := R 2 f ε (t, x, v)dv - R 2 f (t, x, v)dv, x ∈ Ω ρ ε v,2 (t, x) -ρ v,2 (t, x) := R 2 |v| 2 f ε (t, x, v)dv - R 2 |v| 2 f (t, x, v)dv, x ∈ Ω,
in the L ∞ norm in space, for different values of ε, at time t = 1. In Figure 5, these errors are displayed as a function of ε. In both cases, a first order in ε is obtained. These results confirm the strong convergence of Vlasov-Poisson equation (1.1) towards the asymptotic model (1.2).

ρ ε (t = 1, •) -ρ(t = 1, •) L ∞ (left) and ρ ε v (t = 1, •) -ρv(t = 1, •) L ∞ (right) as a function of ε.
Then, we consider the Root Mean Square (RMS) quantity in this example:

RM S(t) := Ω x 2 1 (ρ ε (t, x) -n i )dx, t ≥ 0, (3.2)
which fast oscillates in time with frequency O(ε -2 ) as 0 < ε 1. We consider this quantity using the UA PIC scheme with the initial data U 2nd ± , for ε = 0.025 and ε = 0.0125. On Figure 6, we plot the time history of RM S using the UA PIC scheme under a comparatively large time step ∆t = 0.025, 0.05, 0.1, together with a reference solution obtained with a resolved time step. From the results, we can see that the UA PIC scheme is able to capture the correct RMS quantity at different time even though the time step is too large to resolve the oscillation waves (the ratio ∆t/ε 2 takes values up to 640). Finally, we perform simulations for a larger final time. To validate our simulations, we consider the asymptotic regime ε = 0.005 so that f ε is close to f solution of (1.2). Using an isotropic (in v) initial condition (3.1) implies that the solution f remains isotropic for all time. Then, the asymptotic model (1.2) can be integrated with respect to v to get the guiding-center model (1.3). For this model, the linear theory enables to determine a priori an instability rate for the initially excited modes of the electric potential ( [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF]). This rate can be compared to our numerical results obtained by the UA PIC solver when ε is small. Then, with ε = 0.005, we consider the dynamics of the following quantity

ξ(t) := ln (φ ε 1,1 )(t) , t ≥ 0, (3.3)
where (φ ε 1,1 ) denotes the Fourier coefficient of φ ε of mode (1, 1). Figure 8 shows the time history of ξ(t) under different values of η and k. We observe that when η is small, the linear phase is larger, and for a given k, the numerical rate is in a very good agreement with the analytical one. In Figure 9 (resp. Figure 7), we study the space (resp. velocity) profiles ρ ε (t, x) -n i (resp. χ ε (t, v) := Ω f ε (t, x, v)dx), at different times. The evolution of ρ ε (t, x) -n i is in very good agreement with the results obtained in the literature [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Frénod | Long time behavior of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF]. We can remark in Figure 7 that there is no dynamics in the velocity direction. As mentioned above, this is explained by the isotropic character of the initial condition (3.1).

Example 2

Then, we apply the UA PIC scheme to study the Vlasov-Poisson equation (1.1) with a non-isotropic (in v) initial data:

f 0 (x, v) = 1 4π (1 + sin(x 2 ) + η cos(kx 1 )) e -(v 1 +2) 2 +v 2 2 2 + e -(v 1 -2) 2 +v 2 2 2
.

(3.4)

As discussed previously, since the initial data is no more isotropic and at least when ε << 1, a non trivial dynamics should be observed in the velocity direction. We choose η = 0.05, k = 0.5 and take N p = 409600, N τ = 16. We compute the solution under ε = 0.005 at different times. The dynamics of χ ε (t, v) := Ω f ε (t, x, v)dx is shown in Figure 10. The profiles of ρ ε are similar to the results in (9), so we omit them for brevity. For small times, the two bumps of χ ε start to rotate, at t ≈ 15, the bumps merge and for large times, the velocity distribution becomes isotropic. These effects can not be observed using the guiding-center model (1.3) since the associated initial condition is the integral with respect to v of (3.4).

Long time test

At last but no least, we study the behavior of the UA PIC scheme with the initial condition (3.1) and for a large final time (t f = 100). We fix the numerical parameters as follows: ∆t = 0.1, N x = 128, N y = 64, N τ = 16, N p = 1638400 and consider different values for ε. In Figure 11, the time history of ξ defined by (3.3) is displayed for different values of ε (ε = 1, 0.1, 0.05, 0.01). Different behaviors can be observed (see also [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF]): when ε is small, the Kelvin-Helmoltz instabilities occurs (up to t ≈ 20) whereas when ε is large, we do not get instability anymore. Note that all the numerical parameters are fixed and only the physical parameter ε changes. Furthermore, for large times it turns out that the two-scale t and τ may no longer be separated variables and numerical instabilities may occur. These instabilities have been observed typically for ε = 0.1 in our numerical tests. To avoid this, we propose to "re-initialize" the numerical solution every t 0 (where t 0 has to be chosen, typically t 0 = 20 for this test) as follows. At time t = t 0 , we compute the original unknowns x(t 0 ), v(t 0 ) according to (2.15) then we apply the general procedure where the initial time is now t 0 : two-scale formulation combined with the suitable preparation (2.26) of the initial data at t = t 0 , in which u 0 ± are now computed from x(t 0 ), v(t 0 ) instead of x(0), v(0). This refreshing process is repeated at all times t = kt 0 , k = 1, 2, . . . . To illustrate the influence of this strategy, we plot in Figure 12 the time history of ξ defined by (3.3). The numerical parameters are the same as previously. Without the refreshing process, we observe that numerical instabilities are developed at t ≈ 80 for the UA PIC approach. However, if we apply the refreshing process every t 0 = 20, the solution is smoothed and the numerical instabilities disappear.

Conclusion

In this work, a two-scale Particle-In-Cell solver is introduced for solving the long time behavior of the four dimensional phase space Vlasov-Poisson system under the influence of a strong homogeneous magnetic field. In the asymptotic regime, the solution exhibits oscillations in time of period of order ε 2 where 1/ε is related to the magnitude of the magnetic field.

This numerical scheme is based on a two-scale formulation in which the slow and fast time scale are separated, and on a suitable choice of the initial data for the two-scale system. Then, this solver is proved to be uniformly accurate of second order (in time) with respect to the parameter ε ∈ (0, 1]. Numerical results confirm the good behavior of the scheme in different regimes and considering isotropic and non isotropic initial conditions.

The strategy is able to handle other scaling relevant for strongly magnetized plasmas such as the so-called finite Larmor radius regime. However the case of a non homogeneous magnetic field is still an issue, since oscillations in space and velocity have to be considered.

Figure 1 :

 1 Figure 1: Temporal error of the proposed PIC method in ρ ε with respect to ∆t and ε: results of U 2nd ± (first row); results of U 1st ± (second row); results of U 0th ± (third row).

Figure 2 :

 2 Figure 2: Temporal error of the proposed PIC method under U 2nd ± : results in ρ ε v,1 (first row) and results in ρ ε v,2 (second row) with respect to ∆t and ε.

Figure 3 :

 3 Figure 3: Error of with respect to Nτ under different ε.

Figure 4 :

 4 Figure 4: Time history of the relative error of the total energy for the UA PIC scheme, with ∆t = 0.05, 0.1.

Figure 5 :

 5 Figure 5: Convergence of the Vlasov-Poisson equation (1.1) towards the asymptotic model (1.2):ρ ε (t = 1, •) -ρ(t = 1, •) L ∞ (left) and ρ ε v (t = 1, •) -ρv(t = 1, •) L ∞(right) as a function of ε.

Figure 6 :

 6 Figure 6: Time history of RM S(t) for ε = 0.025 (left) and ε = 0.0125 (right) by using UA PIC with ∆t = 0.025, 0.05, 0.1.

Figure 7 :

 7 Figure 7: Contour plot of quantity χ ε (t, v) at different t with ε = 0.005.

  *0.4*t-5.885

Figure 8 :

 8 Figure 8: Time history of the quantity ξ(t) (3.3) under different η with ε = 0.005, k = 0.5 (first row) and k = 0.4 (second row).

Figure 9 :

 9 Figure 9: Contour plot of quantity ρ ε (t, x) -n i at different t with ε = 0.005.

Figure 10 :Figure 11 :

 1011 Figure 10: Contour plot of quantity χ ε (t, v) at different t with ε = 0.005.

Figure 12 :

 12 Figure 12: Time history the quantity ξ(t) (3.3) for ε = 0.1 without refreshing (left) and with refreshing every t 0 = 20 (right).

Acknowledgements

This work was supported by the French ANR project MOONRISE ANR-14-CE23-0007-01. N. Crouseilles and M. Lemou are supported by the Enabling Research EUROFusion project CfP-WP14-ER-01/IPP-03.