Bastien Le Gloannec
email: bastien.le-gloannec@univ-orleans.fr

Nicolas Ollinger
email: nicolas.ollinger@univ-orleans.fr

Knight tiles: particles and collisions in the realm of 4-way deterministic tilings

Keywords: Deterministic tiles, Domino Problem, Tiling Problem, Expansive Subdynamics, Particles, Wang tiles

Particles and collisions are convenient construction tools to compute inside tilings and enforce complex sets of tilings with simple tilesets. Locally enforceable particles being incompatible with expansivity in the orthogonal direction, a compromise has to be found to combine both notions in a same tileset. This paper introduces knight tiles: a framework to construct 4-way deterministic tilings, that is tilings completely determined by any infinite diagonal of tiles, for which local particles and collisions with many slopes can still be constructed while being expansive in infinitely many directions. The framework is then illustrated by an elegant yet simple construction to mark a diagonal with a 4-way deterministic knight tileset.

A tiling is a coloring of the discrete plane Z 2 assigning a tile from a finite tileset to each position so that the local tiling rules associated to the tileset are satisfied for every tile in its neighborhood. Wang tiles provide a convenient and universal syntactic description of tiling rules: a Wang tile is a unit square with colored edges and the tiling rule requires adjacent squares to share a same color along their common edge -tiling constraints are readable directly on the tiling. Starting with the study of the Domino Problem [START_REF] Berger | The Undecidability of the Domino Problem[END_REF][START_REF] Robinson | Undecidability and nonperiodicity for tilings of the plane[END_REF], computations have been successfully embedded into tilings to enforce computational phenomena and prove undecidability and complexity results. Following what is done in the case of cellular automata, particles and collisions have been successfully used as a tool to transmit information quanta through tilings and mark positions.

A tiling is 4-way deterministic if any infinite diagonal strip of tiles uniquely determines the whole tiling. A syntactic way to enforce such a property on a tileset is to require any consecutive pair of colors along its edges to uniquely determine a tile. Kari and Papasoglu [START_REF] Kari | Deterministic aperiodic tile sets[END_REF] were able to prove that it is still possible to construct aperiodic tilesets under such constraint as are every family of tilings generated by substitution systems [START_REF] Gloannec | Substitutions and strongly deterministic tilesets[END_REF]. Ten years later, Lukkarila [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF] showed that Turing machines can still be simulated in this setting and proved that the Domino Problem is still undecidable for 4-way deterministic tilesets. The construction is subtle and quite involved, building on previous construction by Kari and Papasoglu [START_REF] Kari | Deterministic aperiodic tile sets[END_REF]. The construction has to be quite involved in particular because such a syntactic constraint enforces the expansivity [START_REF] Boyle | Expansive subdynamics[END_REF] of the tilings in every directions but two, prohibiting the use of locally defined particles that are neither horizontal nor vertical. In this paper, inspired by the partitioned cellular automata [START_REF] Morita | Reversible simulation of one-dimensional irreversible cellular automata[END_REF], we propose to replace Wang tiles with knight tiles, keeping a syntactic condition to enforce 4-way determinism and gaining back particles and collisions.

Similar to a Wang tile, a knight tile is a colored unit square sharing colors not only with its 4 direct neighbors but also with the 8 neighbors located at a chess knight move. Determinism along one diagonal is then simply expressed as a syntactic condition on tuples of colors. While keeping expansive in infinitely many directions, knight tiles can be combined to construct particles along slopes that are not too steep. The interest of this family of tilesets is illustrated by considering two problems handled with complicated constructions in Lukkarila [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF].

The first construction is the key ingredient to the proof of the undecidability of the Domino Problem: solving the seeded Domino Problem by simulating a Turing machine.

The second construction is a key ingredient for the previous construction in the classical case: marking a diagonal -a slope that cannot be obtained directly by knight tiles using a particle. Indeed, Lukkarila [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF] asks "Could there be a significantly simpler tile set for drawing a single diagonal line 4-way deterministically?" The construction provided here is significantly simpler, combining time-symmetry and the infinite Firing Squad technique by Kari [START_REF] Kari | Rice's theorem for the limit sets of cellular automata[END_REF] with properties of the Thue-Morse substitution to enforce 4-way determinism.

Preliminary definitions

Wang tiles and deterministic tilesets Given a finite alphabet of colors C, a Wang tile is an oriented unit square tile with one color on each side. Formally, it is a quadruple t ∈ C 4 whose four components are identified with the four directions {w, s, e, n} and, for convenience, denoted as t w , t s , t e and t n . A Wang tileset τ is a finite set of Wang tiles. A tiling by τ is a map T : Z 2 → τ associating a tile of τ to each cell of the discrete plane Z 2 such that two adjacent tiles (for the 4-connectedness) share the same color on their common edge. More formally, a tiling T satisfies the following constraints, for all (x, y)

∈ Z 2 : T (x, y) w = T (x - 1, y) e , T (x, y) s = T (x, y -1) n , T (x, y) e = T (x + 1, y) w , T (x, y) n = T (x, y + 1) s .
The set of all tilings by a tileset τ is denoted as X τ .

The tileset τ is ne-deterministic if for all couples of tiles (w, s) ∈ τ 2 , there exists at most one tile t ∈ τ simultaneously compatible to the west with w and to the south with s: w e = t w and s n = t s . In the case of Wang tiles, one can equivalently say that any tile t ∈ τ is uniquely identified by its couple of colors (t w , t s ). {sw, se, nw}-determinism is defined symmetrically. A tileset is 4-way deterministic if it is simultaneously sw, se, ne and nw-deterministic.

Subshifts and expansiveness Given a finite alphabet Σ, a Σ-coloring (or simply coloring in the absence of ambiguity) of the discrete plane is a map c : Z 2 → Σ. For all u ∈ Z 2 , we define the translation σ u over colorings by

σ u (c)(x) = c(x -u) for all c ∈ Σ Z 2 and x ∈ Z 2 . A coloring c ∈ Σ Z 2 is periodic of period p ∈ Z 2 if σ p (c) = c. The set Σ Z 2 of
all Σ-colorings is endowed with the product (over Z 2 ) topology of the discrete topology (over Σ). A subshift X is a (topologically) closed and translation-invariant (∀u ∈ Z 2 , σ u (X ) = X ) subset of Σ Z 2 . In particular, the set of tilings by a tileset τ is a subshift of τ Z 2 .

For all slope α ∈ R ∪ {∞}, let us denote as l α the real line of slope α going through the origin. For all radius ρ > 0, let us define L α (ρ) = (l α +[-ρ, ρ]u)∩Z 2 for u ∈ R 2 a unit vector orthogonal to l α . L α (ρ) is the discrete thick line of slope α and width 2ρ centered on the origin. A slope α is a direction of expansiveness of a subshift X ⊆ Σ Z 2 if there exists a radius ρ > 0 such that for all x, y ∈ X , x| Lα(ρ) = y| Lα(ρ) =⇒ x = y. In particular, if τ is a 4-way deterministic Wang tileset, X τ is expansive in (at least) all directions of R\{0}.

From Wang tiles to knight tiles

The notion of 4-way determinism introduces a very strong constraint on the tileset. Indeed, it strongly limits its capacity to construct particles locally. To mix particles and determinism, one might loosen the constraints by generalizing the notion of local determinism to a broader radius r ≥ 1 of determinism.

Radius of determinism A tileset τ is ne deterministic with radius r if for all valid (i.e. containing no tiling error along its inner edges) (2r + 1) × (2r + 1) square pattern by τ , the center tile is perfectly determined by the 2r tiles at positions (1, 2r), (2, 2r -1), . . . , (2r -1, 1). Determinism with radius r in the three other diagonal directions (sw, se and nw) is defined symmetrically. The tileset is 4-way deterministic with radius r if it is simultaneously deterministic with radius r in the four diagonal directions.

The different radiuses can be compared through the expansiveness of their tilings as pointed out by the following proposition. Proposition 1. If τ is a 4-way deterministic tileset with radius r, then X τ is (at least) expansive in directions -r r-1 , -r-1 r ∪ r-1 r , r r-1 .

The expansiveness directions are actually tightly related to the particles that can be locally realized. Expansiveness is, up to a change in the radius ρ, independent of the profile chosen for the slope. At a fixed radius of determinism, the ability to mark a line however depends on its profile. Let us denote the horizontal and vertical unit vectors of Z 2 as e 1 = (1, 0) and e 2 = (0, 1). For any slope α ∈ Q ∪ {∞}, a (periodic) profile for α is a finite sequence (u 0 , . . . , u k-1 ) ∈ {e 1 , e 2 } k such that α = i ui•e2 i ui•e1 . Given a profile P = (u 0 , . . . , u k-1 ) for α, we define the associated subshift L P ⊆ {0, 1} Z 2 as the subshift generated by (i.e. the smallest subshift containing) the configuration c P containing a discrete 4-connected line of slope α drawn following the profile P , i.e. formally defined by c P (x) = 1 if x = 0≤i≤n u i mod k or x = --n≤i<0 u i mod k for some n ≥ 0, and c P (x) = 0 for all other x ∈ Z 2 . Note that L P exactly contains all translated versions of c P , plus the blank coloring (all 0s) for compactness reasons. A coloring c ∈ Σ Z 2 , where Σ contains an identified blank color c 0 ∈ Σ, is a particle of direction α ∈ Q ∪ {∞} if it is periodic in direction α (i.e. c admits a periodicity vector (p x , p y ) ∈ Z 2 such that α = py px ) and ultimately constant1 equal to c 0 in any other direction. This equivalently means that all non-blank cells of c are contained in L α (ρ) + ∆ for some ρ > 0 and ∆ ∈ Z 2 . We say that a slope α ∈ Q ∪ {∞} is locally realized by a tileset τ if there exists a projection π : τ → {0, 1} (naturally extended to colorings π : τ Z 2 → {0, 1} Z 2 ) such that there exists a profile P for the slope α such that π(X τ ) = L P and every tiling of X τ is a particle. The previous definition introduces a reinforced notion of soficity for L P .

The following results conclude our remarks on radiuses of determinism.

Proposition 2. If a slope α is locally realizable by a tileset τ then α is not a direction of expansiveness for X τ .

Proposition 3. The slopes that are locally realizable by 4-way deterministic tilesets at radius r are exactly

Q ∩ -∞, -r r-1 ∪ -r-1 r , r-1 r ∪ r r-1 , +∞ .
Knight tiles Whereas classical 4-way determinism is a purely syntactic property of the tileset, to check determinism at radius r one has to consider all tilings of (2r + 1 × 2r + 1) squares. We introduce deterministic knight tiles as a convenient and purely syntactic notion of determinism at radius 2 that can be checked directly on the tiles. Given a finite alphabet C of colors, a knight tile is formally a 12-tuple of C 12 and a knight tileset is a finite set of knight tiles. For convenience, each of the twelve components of a knight tile will be identified by a direction among {w, s, e, n, ws, sw, se, es, en, ne, nw, wn} and for a knight tile T ∈ C 12 , we will for instance denote by T w ∈ C the corresponding w color. A tiling of the discrete plane by a knight tileset τ is a map t : Z 2 → τ satisfying the following constraints, for all (x, y) ∈ Z 2 : t(x, y) andt(x, y) wn = t(x -2, y + 1) es . This means that each tile of a tiling shares each of its different color components with one of the twelve neighboring tiles represented on the figure 1(a).

w = t(x -1, y) e , t(x, y) s = t(x, y -1) n , t(x, y) e = t(x + 1, y) w , t(x, y) n = t(x, y + 1) s , t(x, y) ws = t(x -2, y -1) en , t(x, y) sw = t(x -1, y -2) ne , t(x, y) se = t(x + 1, y -2) nw , t(x, y) es = t(x + 2, y -1) wn , t(x, y) en = t(x+2, y+1) ws , t(x, y) ne = t(x+1, y+2) sw , t(x, y) nw = t(x-1, y+2) se ,
A knight tileset τ is ne-deterministic if any tile t ∈ τ is uniquely identified by its quadruple of colors (t wn , t w , t s , t se ), i.e. there is at most one tile in the set that is compatible with the four colors pointed out on the figure 1(b). Determinism A knight tileset is 4-way deterministic if it is simultaneously deterministic in the four diagonal directions sw, se, ne and nw.

The rest of the paper is dedicated to advocate for the use of this notion of radius-2 determinism by showing that constructions that are painful to handle in the classical 4-way deterministic setting can be treated with particles and collisions in the 4-way deterministic knight setting.

The seeded Knight Domino Problem is undecidable

Simulating a Turing machine with classical 4-way tiles requires a complicated machinery [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF] involving particular aperiodic tilesets. Knight colors provide enough flexibility to handle it directly and prove the undecidability of the Domino Problem with a seed tile in a classical way.

Problem 1 (Domino Problem with a seed tile). Given a tileset τ and a specified seed tile t 0 ∈ τ , does there exist a tiling of Z 2 by τ using (at least once) the tile t 0 ?

The reader is assumed familiar with Turing machines. As every Turing computation can be made reversible [START_REF] Lecerf | Machines de Turing réversibles[END_REF][START_REF] Bennett | Logical reversibility of computation[END_REF], we directly work with reversible Turing machines. A reversible Turing machine (RTM) is a 5-tuple (Σ, Q, ← → q i , F, δ) where Σ is the tape (finite) alphabet, Q the finite set of states. Before defining the remaining elements of the tuple, let us state that the head always moves at each transition and we denote its two possible moves by {←, →}. Let us also define ← → Q = Q × {←, →} whose elements will be written -→ q (resp. ←q ) to denote (q, →) (resp. (q, ←)). Then we define the partial injective transition map δ :

← → Q × Σ → ← → Q × Σ. Finally ← → q i ∈ ← → Q is the initial (oriented) state and F ⊂ ← → Q is the set of final (oriented)
states (and we will assume that δ is not defined on states of F ). This is an acceptable model for reversible Turing machines that will be furthermore well-fitted to our purpose.

The tileset described on figure 2, where colors are conveniently represented by arrows that must be go uninterrupted across tile edges, is a 4-way deterministic knight tileset that simulates a given RTM (Σ, Q, ← → q i , F, δ). Firstly note that the only knight color components used in this simulation are {sw, se, ne, nw} hence there is no ambiguity on the knight colors represented on the figure 2.

To be able to use efficiently these knight colors for determinism, one must slow down the simulation. To that purpose, each transition of the Turing machine is decomposed into three steps in this fixed order: 1. a state transition step where the transition is done but the head does not yet move, though its move is already contained in the state of ← → Q (the corresponding tiles are represented on the third line of the figure); 2. a waiting step where nothing happens (second line of the figure); 3. a move step where the head finally moves (first line of the figure). The Wang constraints to enforce this order for the steps of the computation process are not explicitely represented on the tiles of figure 2, but it is a simple counting modulo 3 in the vertical colors. Assuming that a Turing computation is correctly initialized in a tiling, then each line of the space-time diagram of the Turing machine can be read once every three lines of the tiling.

For all a, b ∈ Σ, ← → q ∈ ← → Q : a a ← -q -→ q b b ← -q a ← -q a ← -q ← -q a -→ q a -→ q -→ q -→ q ← -q b b -→ q a a ← → q a ← → q a ← -p δ -1 ( ← → q , a) = ( -→ p , b) ← → q a ← → q a -→ p δ -1 ( ← → q , a) = ( ← -p , b) a a ← → q a ← -r b ← -r δ( ← → q , a) = ( ← -r , b) ← → q a -→ r b -→ r δ( ← → q , a) = ( -→ r , b)
Fig. 2. Deterministic knight tiles simulating a RTM One can enforce the initialization of the Turing computation using a seed tile and the tiles of the figure 3(a) that force a blank tape containing a unique Turing machine head to appear in the tilings. The labels "≤ -2", "-1", "0", "1" and "≥ 2" appearing on the tiles denote some specifications for the horizontal colors: the tile with index 0 is only compatible to the left (resp. right) with the tile with index -1 (resp. 1) and, to the left (resp. right) of the -1 (resp. 1) tile, only the "≤ -2" (resp. "≥ 2") tile can appear. However, these tiles cannot directly be added to the previous simulation tile set without losing the 4-way determinism. Indeed, it not difficult to see that, as is, predicting the tiles with index -1 or 1 is not possible in every direction. To solve that problem, one can add a very simple layer of information represented by the Wang tiles of the figure 3(b): a blank tile g 0 and a tile g 1 carrying a vertical "ghost" signal (that must be vertically preserved along columns of the tilings). Denoting as τ T the Turing simulation tiles of the figure 2,τ 1 i the three initialization tiles with indices {-1, 0, 1} of the figure 3(a),τ 2 i the two initialization tiles with labels "≤ -2" and "≥ 2" of the figure 3(a) and b the blank tile of the figure 3(a), we define the two-layered knight tileset τ = ((τ T ∪ {b}) × {g 0 , g 1 }) ∪ (τ 1 i × {g 1 }) ∪ (τ 2 i × {g 0 }), which can be interpreted as a standard knight tileset on couples of colors of each layer. In a tiling, the three consecutive columns containing the initialization tiles with indices {-1, 0, 1} are the only one carrying a ghost signal (tile g 1 ) on the second layer. This allows to predict these tiles in every direction, hence the knight tileset τ is 4-way deterministic. As the transition tiles are not defined on final states, τ tiles the plane if and only if the simulated Turing machine does not halt on the blank tape. It is then straightforward to derive the following result.

Theorem 1. The Domino Problem with a seed tile is undecidable for 4-way deterministic knight tilesets. Considering deterministic knight tiles allows to easily build particules that are not realizable in the usual deterministic Wang framework (where the only realizable slopes are 0 and ∞ by proposition 3) as illustrated by the following result.

B ≤ -2 B -1 -→ qi if ← → qi = ← - qi or B -1 if ← → qi = -→ qi ← → qi B 0 seed B 1 if ← → qi = ← - qi B 1 ← - qi or if ← → qi = -→ qi B ≥ 2
Proposition 4. The slopes that are locally realizable using 4-way deterministic knight tiles are exactly

Q ∩ [-∞, -2] ∪ -1 2 , 1 2 ∪ [2, +∞] .
Our main construction will rely on these realizable slopes to mark a diagonal in a lighter way than the construction of [START_REF] Lukkarila | The 4-way deterministic tiling problem is undecidable[END_REF] in the Wang case. The general idea is to use a fundamentaly one-dimensional hierarchical structure of signals similar to an infinite version of Minsky's classical solution to the Firing Squad synchronization problem on cellular automata, as what is done for instance in [START_REF] Kari | Rice's theorem for the limit sets of cellular automata[END_REF].

This structure has moreover to be reversible so that the tileset can be made deterministic in two opposite directions, say ne and sw. The structure also has to be compatible with determinism for knight tiles in the two orthogonal directions nw and se. The reader might keep in mind that only knight colors {ws, sw, en, ne} will be used. Thus the tileset will be deterministic at radius 1 (classical Wang tiles determinism) in directions ne and sw, while it will be deterministic at radius 2 (using knights) in directions nw and se.

The general structure is described on figure 4 where the red dots represent the diagonal to be marked, referred to as the fire line. For convenience, the four kinds of signals are named H, H , V and V according to the figure. Their respective slopes are 0, 1 3 , 2 3 and ∞. A binary hierarchical structure is used to mark some points along successive front lines until the granularity of space allows all the marked points of a front to be sufficiently close to decide by a local rule to fire at that point (and consequently mark the diagonal). Let us enumerate the front lines by their rank starting at 0 for the fire line. The points of the front of rank k are regularly arranged and two consecutive points are separated by 2 k -1 non-marked positions. The marked points of the front of rank k will be referred to as pillars, which correspond to the points that are marked on the front of rank k + 1 (i.e. points that were already marked seen things as a firing squad with time running towards ne), and middles, which exactly correspond to the middle positions between pillars. Referring to the figure 4 for the orientation, building a local rule for the construction to be deterministic in the ne direction is not particularly difficult and goes straightforward using the represented signals. By contrast, making this construction deterministic in the sw direction is more challenging as one has to be able to locally and deterministically distinguish between pillars and middles. For this, one need an infinite source of well-structured alternating bits sequences. A way to get such information properly is to resort to a sequence that is a fixed-point of a well-chosen substitution.

Let us define the Thue-Morse substitution (see [START_REF] Fogg | Substitutions in Dynamics, Arithmetics and Combinatorics[END_REF]) s : {0, 1} → {0, 1} * by s(0) = 01 and (1) = 10. We naturally extend the definition of s to finite words of {0, 1} * (resp. infinite words of {0, 1} N ) defining, for all u ∈ {0, 1} * (resp. u ∈ {0, 1} N ), s(u) 2i+k = s(u i ) k for all i ∈ {0, . . . , |u| -1} (resp. i ∈ N) and k ∈ {0, 1}. The Thue-Morse word T ∈ {0, 1} N is the unique fixed-point of s starting with the letter 0: T = lim n→+∞ s n (0) that exists as s(0) 1 = 0 and is infinite as |s(0)| > 1. As the Thue-Morse word is a fixed-point for s, it can be unsubstituted into itself by s placing bars every two letters starting with a bar at position 0. A given factor v of the Thue-Morse word is said even if there exists a position p ∈ N even such that T p • • • T |v| = v, which means that there is a bar just before its first letter v 0 (i.e. between T p-1 and T p in the Thue-Morse word) in one of the possible decompositions of v. Complementarily, a factor is said odd if it appears at an odd position p in the Thus-Morse word, which means that the bar is between v 1 and v 2 in a decomposition of v. A factor can simultaneously be even and odd if it appears in the Thue-Morse word at both even and odd positions: for instance 011 appears at positions 1 and 7. However,

H V V H fi r e
li n e f r o n t 0 f r o n t 1 f r o n t 2 Fig. 4. General scheme of the structure the following lemma points out the fact that a sufficiently large factor cannot be simultaneously even and odd. Lemma 1. Any factor of T of length at least four is either even or odd, which means that it admits a unique decomposition.

Hence it is sufficient to look at factors of size four to be able to determine unambiguously their alignment. Observe that moving letter by letter a window of size four over the Thue-Morse word, one sees an alternating sequence of even/odd factors. This will constitute the alternating bits sequence we require. Now that all the required objects are defined, each signal of the construction will henceforth carry a factor of size four of the Thue-Morse sequence. The signals H and V only carry the even words while H and V only carry odd words. When two signals meet on a marked point, we require them to carry the same factor u (|u| = 4) and we substitute it by s (|s(u)| = 8) to derive the four factors carried by the four "outgoing" (seeing time going ne) signals: the produced type V signal must carry s(u)

1 • • • s(u) 4 , the type V must carry s(u) 2 • • • s(u) 5 , the type H must carry s(u) 4 • • • s(u) 7 and the type H must carry s(u) 5 • • • s(u) 8 .
Firstly remark that knowing s(u), one uniquely deduces u. Although there is no "middle signal" to carry s(u) 3 • • • s(u) 6 , it can easily be deduced from the four other derived words we dispose of.

This defines a coherent structure. Indeed, let us consider any bi-infinite word of {0, 1} Z such that every factor appears in the Thue-Morse sequence2 , "written" along a front line of rank k > 1 in the following sense: two consecutive

0 0 1 0 1 1 0 0 1 1 1 1 0 0 1 1 0 1 0 0 1 0 0 1 0 0 1 1 0 0 1 0 0 1 1 1 0 1 1 0 1 1 0 0 1 1 1 0 0 1 0 1 0 1 1 0 1 1 0 0 0 0 1 0 1 1 0 1 1 0 1 1 0 1 0 0 1 1 0 0 1 0 0 1 0 0 1 0 1 1 0 0 1 0 1 0 1 0 0 1 0 0 1 0 1 0 0 1 1 0 1 0 1 0 0 1 0 0 1 0 1 0 0 1 0 1 1 0 1 0 0 1 0 0 1 0 1 1 0 1 1 0 0 1 1 0 1 1 0 1 1 0 1 0 0 1 1 0 0 1 1 0 1 1 0 1 1 0 1 0 0 1 0 0 1 0 1 0 1 1 0 1 1 0 1 0 1 1 0 0 0 0 1 0 1 1 0 1 1 0 1 0 1 0 0 1 0 0 1 0 0 1 0 1 1 0 0 1 1 0 1 1 0 0 1 0 0 1 0 1 1 1 0 1 0 0 1 1 0 0 1 0 0 1 0 1 0 0 1 1 1 0 1 0 0 0 0 1 1 0 1 1 0 0 1 0 1 1 0 1 1 0 0 1 1 0 1 1 0 0 0 1 0 0 1 0 1 0 0 1 0 0 1 1 0 1 0 1 0 0 1 0 0 1 1 1 0 1 0 0 1 1 0 1 0 1 1 0 1 0 0 1 0 0 1 0 1 1 0 1 1 0 1 0 0 1 1 0 1 0 1 1 0 1 0 0 1 1 0 1 0 1 1 0 1 0 1 1 0 1
Fig. 5. Marking a diagonal marked points hold two consecutive factors of size 4 with a one-letter shift. The previously described local mecanisms of substitution/unsubstitution on marked points enforce the same word (as it is a fixed-point) to appear on all the other front lines and to be used as a ressource to properly alternate between middles and pillars on a front line, coherently with other lines so that a proper global binary hierarchical structure is enforced.

A realization of this structure using knight tiles is given on the figure 5 where the colors are represented using some convenient signal notations. The knight colors are not represented on this figure but are all illustrated on the figure 6. H and V signals, of respective slopes 0 and ∞, carry even words and meet on pillars (solid black squares on the figures). The pilar tiles force their two incoming H and V signals to carry the same factor u. As the four outcoming factors are derived from this common word, a pillar tile is perfectly determined by u, which is explicitely written next to each pillar on the figure 5. Similarly, H and V signals, of respective slopes 1 3 and 2 3 realized using particular profiles for reasons that will be exposed later and two different types of knight colors to distinguish the two different types of angles that appear in that profile as illustrated by figures 6(a) and 6(b), carry odd factors. These signals meet on middles, that are represented by gray squares on the figures. The middle tiles force their two incoming H and V signals to carry the same factor u. As the four outcoming factors are derived from this common word, a middle tile is perfectly determined by u, which is explicitely written next to each middle on the figure 5. We do not dispose of enough space to formally describe the whole tileset. Although most of it is a rather straightforward translation of the previously exposed scheme (figure 4 and its symmetric along the fire line) into a tileset, several points certainly require explanations. This is what the remaining paragraphs are dedicated to.

(a) H signals (b) V signals (c) Pillar (d) Middle (e) Fire
First, when a H signal meets a V , they do not simply cross each other, there is actually a collision occuring (grey dot on the figures) and there are actually two kinds of H/V signals. When two signals of the first kind meet, they cross each other and by the way change into the second kind (darker color on the figures). When two signals of the second kind meet, they produce a pillar that will in turn produce two signals of the first kind. This collision has to be treated particularly when it is the last collision before the fire line, which can be detected as the collision tile also has to contain some H /V signals when it is near the fire line (see figure 6(e)). In that case, the first kind of H/V signals is changed into a third kind that produces a fire tile (instead of a pillar) when meeting.

The case of H and V signals requires some careful attention as we need to use knight colors for them to be deterministic in the nw and se directions. The knight colors used basically contain the type of signal and the word it carries. However a knight must absolutely not cross a front line as in this case the determinism in one of the orthogonal ne or sw directions would fail. Without entering into the details, our solution consists in using the (e 1 , e 1 , e 2 , e 1 , e 1 , e 1 , e 1 , e 2 ) profile (alternating between 2 and 4 horizontal steps) such that one can make sure that the tiles that appear at both end of a signal are tiles with blank knight colors (hence no knight crosses the front line). We use the symmetrical profile for V signals. All situations are illustrated by figure 6.

The fire line must appear when marked points are sufficiently close (and as for the front line, no knight should cross the fire line). We resort for that to several ad-hoc tiles that appear only near the fire line, as this was already illustrated by the third kind of H/V signals. This is illustrated by the figure 6(e).

On the other side of the fire line, say the ne side, the hierarchical structure is symmetrically dismantled. This is done by duplicating all the colors used in the sw part and defining a ne tile from each sw tile by applying a symmetry along the nw-sw diagonal on the tile and replacing each sw color by its ne duplicate. Only the fire tiles, that make the junction between boths sides, simultaneously hold sw and ne colors. That way, the resulting tileset is time-symmetric: it is its own "inverse" up to a swap in colors.

To conclude, let us denote as κ the previously described knight tileset. The following theorem underlines the properties κ was builded for.

Theorem 2. κ is a 4-way deterministic knight tileset and its set of tilings projects onto the diagonal subshift, i.e. {0, 1}-colorings containing at most one diagonal of 1 and 0 everywhere else.

Appendix: Proofs and complementary material

The complementary material for the reviewers follows the order of the article.

Proof of proposition 1 A perfectly formalized proof of this result would require to define some formal tools to properly manipulate the objects we are dealing with, like the coding relation of [START_REF] Boyle | Expansive subdynamics[END_REF] (definition 3.1). For brevity and simplicity, we will only give a slightly informal yet convincing presentation of our arguments.

Proof sketch. First observe that, for α ∈ R∪{∞}, if there exist ρ, ε > 0 such that L α (ρ) L α (ρ+ε) and for all x ∈ X τ , x| Lα(ρ) perfectly determines x| Lα(ρ+ε) , then X τ is expansive in direction α (as, "iterating" the argument, x| Lα(ρ) perfectly determines x). One can refer to lemma 3.3 of [START_REF] Boyle | Expansive subdynamics[END_REF] for more details.

For a tileset se and ne-deterministic at radius r, the two red cells of the subfigures of the figure 7 are uniquely determined by the 2r white ones. For any slope α ∈ r-1 r , r r-1 , i.e. strictly between the two limit slopes in blue on the figures, there exists an appropriate ρ 0 > 0 such that, up to a translation, the red cells (i.e. their centers) can be made outside the "expansiveness window" L α (ρ 0 ) while the (centers of the) white ones all are inside. This argument applied along any expansiveness window of radius ρ ≥ ρ 0 ensures that it can be locally extended on both sides using the local rule of determinism. Said otherwise, as illustrated by the figure 8 where the rule of determinism is represented by the red lozenge lying on the line of slope 1, there exists ε > 0 such that L α (ρ) L α (ρ+ε) and for all x ∈ X τ , x| Lα(ρ) perfectly determines x| Lα(ρ+ε) which implies the expansiveness in direction α.

The other interval of directions of expansiveness is obtained symmetrically considering the two other directions of determinism.

Proof of proposition 2

Proof. For any ρ > 0, Consider the tiling that is projected to the blank coloring (all 0s) and a tiling that projects onto a line of slope α sufficiently far from the origin so that only blank tiles are visible in the "expansiveness window" L α (ρ). Both tilings coincide on L α (ρ) but are distinct, hence X τ cannot be expansive in direction α.

Proof of proposition 3

Proof. The realizability part is proved in details for radius 2 in the more constraining case of knight tiles (keep the same profile, simply remove the additional knight colors that are not required anymore as the local rule of determinism can already "see" the cells "emitting" these colors) in the proof of proposition 4. This generalizes straightforwardly to the case of deterministic Wang tiles of arbitrary radius.

The "exactly" part is a consequence of propositions 1 and 2. More detailed explanations for H and V signals The case of H and V signals requires some careful attention. Indeed we need to use knight colors for them to be deterministic in the nw and se directions. The knight colors used basically contain the type of signal and the word it carries. However there is a constraint on the possible positions for a knight color: a knight must absolutely not cross a front line as in this case the determinism in one of the orthogonal ne or sw directions would fail. Indeed, seeing time going ne, this could happen in the two following situations:

1. The knight links a tile from an incoming H or V signal to a tile that is on the ne side of the front line. Then the factor carried by this knight would not be predictable in the sw direction as the four derived factor would not all be "readable" by the deterministic rule. 2. The knight links an outcoming H or V tile to a tile that is on the sw side of the front line. Then this tile would not be predictable in the ne direction.

Using the quite natural (e 1 , e 1 , e 1 , e 2 ) profile (3 horizontal steps for each vertical one) for the slope 1 3 of H signals would make us fall into one of these situations. The adopted solution consists in using the (e 1 , e 1 , e 2 , e 1 , e 1 , e 1 , e 1 , e 2 ) profile (alternating between 2 and 4 horizontal steps). That way, there are exactly three tiles in the signal that do not hold a knight color: the tile that is in the middle of the 3 tiles corresponding to the 2 horizontal steps and the two tiles surrounding the middle of the 5 tiles corresponding to the 4 horizontal steps. We use the symmetrical profile for V signals. Boths cases are illustrated on the figures 6(a) and 6(b). Moreover, by choosing the right starting tile for the emission of a signal on a pillar/middle, one can make sure that when two H and V signals meet to produce a middle, they had always just used that knight-less tile, i.e. the western and southern neighbors of a middle are always the two knight-less tiles of the H and V signal respectively. Said otherwise, for H signals (the V case is symmetric), the marked points always appear at the position of the middle tile in the 5 tiles corresponding to the 4 horizontal steps, hence they are directly surrounded in the signal by tiles that hold no knight color. We prevent that way the first failure as illustrated by the figure 6(d): the knights from the incoming signals do not cross the front line. The choice done also prevents the second failure as illustrated by the figures 6(c) and 6(d): the knight from the outcoming signals do not cross the front line.

Complete description of the tileset The tileset κ contains 1151 tiles that can be classified into the 41 categories of figure 11. The colors are denoted by lines drawn on the tiles that must continue uninterrupted across tile edges, except for the H and V signals where, due to their particular profiles, some discontinuities appear along edges marked by the special colors # 1 and # 2 . The H and V signals can carry any of the 4 even subwords {0101, 1001, 1010, 0110}. The H and V signals can carry any of the 6 odd subwords {1011, 0010, 1100, 0011, 0100, 1101}. The knight colors only appear on tiles containing a H or V signal and carry the same odd subword as the associated H or V signal. On tiles containing several H and V signals, the direction of each knight color easily permits to identify unambiguously the associated H or V signal. Pillar and middle tiles were already comprehensively described in the article (the two incoming signals carry the same subword, the four outcoming words are derived from the substitution of this subword).

Only the tiles of the sw half plane delimited by the fire line are represented on the figure 11. The ne tiles appearing on the other side of the fire line are obtained from the sw tiles by applying the transformation represented on the figure 12 (color swap and symmetry along the nw-sw diagonal of the tiles), where, for each sw color c, we denote by c its duplicate for the ne tileset. The same notation is used on figure 11 for the fire line tiles, which must themselves respect this symmetry as they make the junction between both symmetric tilesets. The transformation is applied to all the tiles to obtain their symmetrical counterpart, except for the fire line tiles (that are fixed points for this transformation).

A large version of the figure 5 is provided on figure 13.

Proof of theorem 2

Proof sketch. The first assertion is a straightforward syntactic verification on tiles. Note that κ is even deterministic at radius 1 (standard Wang tiles determinism) in directions ne and sw(as the knight colors for these directions are not used here).

On the other hand clear, κ tiles the plane (blank tiling) and there exists π : κ → {0, 1} (naturally extended to colorings π : κ Z 2 → {0, 1} Z 2 ) such that π(X κ ) = L -1 , where L -1 denote the subshift generated by the {0, 1}-coloring c such that, for all (x, y) ∈ Z 2 , c(x, y) = 1 if x + y = 0 and c(x, y) = 0 otherwise. Indeed, picking an arbitrarily long factor of the Thue-Morse sequence and writing it with fire tiles on the diagonal of an empty square, as on the figure 5, one easily sees that the local ne and sw determinisms can be applied to turn the whole square into a valid square pattern. Hence there exist arbitrarily large valid square patterns marking their finite diagonal. For compactness reasons, there must exist a tiling of the whole plane containing a marked bi-infinite diagonal. This implies in turn the left-right inclusion by translation invariance. To prove the converse, observe that a tiling containing a fire tile must contain a full fire line and the whole associated structure, hence all "degenerated" tilings are projected to the blank coloring.

Finally, observe that the set of tilings by κ is at least expansive in directions ]-∞, 0[ ∪ 1 2 , 2 .

Blank tile:

H and V signals and collision:

H signal: #1 #1 #2 #2 V signal: #1 #1 #2 #2
V /H and H/V crossings: Pillar, signal emission tiles and middle: Particular collisions and crossings near the fire line:

Fig. 1 .

 1 Fig. 1. Knight tiles and determinism

Fig. 3 .

 3 Fig. 3. Additional knight tiles to initialize the Turing computation

Fig. 6 .

 6 Fig. 6. Knight colors of the tileset

Fig. 7 .Fig. 8 . 1 Fig. 10 .

 78110 Fig. 7. Limit expansive directions for determinism radius r

Fig. 11 .Fig. 12 .

 1112 Fig. 11. 41 types of tiles in κ

Note that all our results would remain true replacing "ultimately constant" by "ultimately periodic" in the definition of a particle.

For topological reasons, such a word exists. The set of all these words is actually a non-empty subshift of {0, 1} Z . All of them are fixed-point for the Thue-Morse substitution.

Example of tiling simulating a Turing computation

Proof of proposition 4

Proof. Firstly recall that the "exactly" part is due to propositions 1 and 2. Also remark that the cases of slopes 0 and ∞ are trivial, the case of the interval -1 2 , 0 is equivalent to the case of the interval 0, 1 2 up to a horizontal symmetry on tiles and the case of the interval [2, +∞[ (resp. ]-∞, 2]) is equivalent to the case of the interval 0, 1 2 (resp. -1 2 , 0 ) up to a symmetry along the sw-ne (resp. nw-se) diagonal on tiles. Hence it is enough to consider the case of a slope α ∈ Q∩ 0, 1 2 . For such a slope, the determinism in directions ne and sw can be obtained without effort (as one directly sees the tiles of the slope in these directions), hence only knight colors of types {ws, sw, en, ne} could be useful here, and we will actually only use {ws, en}. The difficulty is to be able to predict the appearance of the angles of the slope in the nw and se directions. All we need is to provide a profile (considered cyclically) for α whose angles (i.e. transitions e 1 /e 2 ) can be predicted at radius 2 using knight colors. It it not difficult to see that as long as no sub-sequence of the form (e 1 , e 2 , e 1 , e 2 , e 1 , e 2 ) (resp. (e 2 , e 1 , e 2 , e 1 , e 2 , e 1 )) appears in the profile, it is always possible to "see" the line with knight colors of types {sw, en} (resp. {ws, ne}) when going towards nw (resp. se). If we limit ourselves to the single knight color en (resp. ws), we have the more restrictive (e 1 , e 2 , e 1 , e 2 ) (resp. (e 2 , e 1 , e 2 , e 1 )) forbidden pattern. Let us write α = p q with q ≥ 2p > 0 and consider the profile P = (e 1 , e 1 , e 2 , . . . , e 1 , e 1 , e 2 p times , e 1 , . . . , e 1 q-2p times ).

It satisfies the required property and using a special knight color to mark each of the angles of the line, one straightforwardly builds 4-way deterministic knight tiles to locally realize the line.

Parabola We give here without entering into the details an example of a nontrivial system of particles and collisions: the parabola of the figure 10. Such a system would be technically very complex to realize using 4-way deterministic Wang tiles. It is however rather simple to achieve in the 4-way deterministic knight setting. Only the knight colors {sw, se, ne, nw} are used. The tiles can directly be read on the figure.

Proof of lemma 1

Proof. Let us first state one classical fact about the Thue-Morse word: T does not contain any triple factor www where w ∈ {0, 1} * . A proof (and more about the Thue-Morse substitution) can for instance be found in [START_REF] Fogg | Substitutions in Dynamics, Arithmetics and Combinatorics[END_REF] (chapter 5). Also observe that if 00 or 11 appears, then the bar is necessarily in the middle (0|0 or 1|1) as they do not belong to s({0, 1}) = {01, 10} (and this observation is not ambiguous as the triple factors 000 and 111 do not appear). Among all 16 words of {0, 1} 4 , only 10 do not contain a triple factor (and they all actually appear in the Thue-Morse word). All of them contain 00 or 11 (thus are unambiguous) except 0101 and 1010. As they are symmetric, we only consider 0101. If it were odd, i.e. 0|10|1, then we would necessarily have 10|10|10 considering the previous and the next letter of the factor in T , which is contradictory as this gives


S Ébastien Guisset

S T Éphane Brull
email: stephane.brull@math.u-bordeaux1.fr

Emmanuel D'humi Ères
email: dhumieres@celia.u-bordeaux1.fr

Bruno Dubroca
email: dubroca@celia.u-bordeaux1.fr

ASYMPTOTIC-PRESERVING SCHEME FOR THE FOKKER-PLANCK-LANDAU-MAXWELL SYSTEM IN THE QUASI-NEUTRAL REGIME

Keywords: Asymptotic-Preserving scheme, Fokker-Planck-Landau equation, Maxwell equations, quasi-neutral limit, angular M1 moments model

This work deals with the numerical resolution of the Fokker-Planck-Maxwell system in the quasi-neutral regime. In this regime the stiffness of the stability constraints of classic schemes causes huge calculation times. That is why, we introduce a new stable numerical scheme consistent with the transitional and limit models. Such schemes are called Asymptotic-Preserving (AP) schemes in literature. This new scheme is able to handle the quasi-neutrality limit regime without any restrictions on time and space steps. This approach can be easily applied to angular moment models by using a moments extraction. Finally, two physically relevant numerical test cases are presented for the Asymptotic-Preserving scheme in different regimes. The first one shows the efficiency of the Asymptotic-Preserving scheme in the quasi-neutral regime whereas the second one on the contrary corresponds to a regime where electromagnetic effects are predominant.

1. Introduction. This work deals with non-homogeneous collisional plasmas described by a kinetic model. The plasma is considered as a mixture of electrons and ions. Each species is characterised by its distribution function which corresponds to particles density in the phase space. In this work the Fokker-Planck-Landau equation describes electron transport and consider binary collisions between particles whereas Maxwell equations are used to describe the evolution of electromagnetic fields. For the sake of simplicity, we assume that the plasma consists of electrons and one ion species considered as immobile. This approximation is relevant due to the important mass of ions compared to the electrons mass. This means our model is valid on time scales during which the ions motion can be neglected.

For the study of collisional processes, the two important physical scales are the mean free path and the electron-ion collision frequency. The mean free path represents the average distance travelled by an electron between two collisions with an ion. The electron-ion collision frequency represents the number of electron-ion collision per unit of time. When the electronic plasma period is very small compared to the electron-ion collisional period and the Debye length is very small compared to the mean free path, the plasma is designated as quasi-neutral and the Maxwell-Gauss (also called Maxwell-Poisson) and Maxwell-Ampere equations degenerate into algebric equations on collisional time scales.

Therefore to handle this type of situation a new class of methods, called Asymptotic-Preserving (AP) methods has been developed. These methods have been introduced firstly by Shi Jin ([23]) in the context of diffusive limits for kinetic equations. Consider a system (S α ) depending on a parameter α and (S 0 ) the corresponding limit system when α tends to zero. In our case α is the ratio between the Debye length and the mean free path. A numerical scheme with time step ∆t and space step ∆x is called Asymptotic-Preserving in the limit α tends to zero for the system (S α ) if the scheme is stable independently of the values taken by α and if the limit scheme obtained for α = 0 is consistent with the limit problem (S 0 ). In this work the system (S α ) corresponds to the Fokker-Planck-Landau-Maxwell system and (S 0 ) corresponds to the Fokker-Planck-Landau-Maxwell system in the quasi-neutral limit. This regime has been already studied in the context of fluid models ( [START_REF] Crispel | A plasma expansion model based on the full Euler-Poisson system[END_REF][START_REF] Degond | Numerical approximation of the Euler-Poisson-Boltzmann model in the quasineutral limit[END_REF]). In ( [START_REF] Crispel | An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasi-neutral limit[END_REF]), the authors considered a two fluid isentropic Euler system coupled with the Poisson equation. It is shown that the Maxwell-Poisson equation can be reformulated into an elliptic equation which does not degenerate at the quasi-neutral limit. In ( [START_REF] Crispel | An Asymptotic Preserving scheme for the Euler equations in a strong magnetic field[END_REF]), this approach is generalised to the Euler-Maxwell model with a strong magnetic field. A kinetic model consisting in a two fluid Vlasov-Poisson system has also been investigated in ([12]). In ( [START_REF] Degond | Numerical approximation of the Euler-Maxwell model in the quasineutral limit[END_REF]), an Asymptotic-Preserving scheme is proposed for the Euler-Maxwell system in the quasi-neutral regime. The Maxwell equations are reformulated to enable the computation of the electrostatic field even in the limit regime. The development followed to express the electric field is well known in physics ( [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF], [START_REF] Braginskii | Reviews of plasmaphysics[END_REF]).

The present paper deals with the construction of an Asymptotic-Preserving scheme for the Fokker-Planck-Landau-Maxwell system in the quasi-neutral limit. The strategy adopted is similar to the one in ( [START_REF] Degond | Numerical approximation of the Euler-Maxwell model in the quasineutral limit[END_REF]), nevertheless to our knowledge, it is the first time that such schemes are considered for kinetic models with true collision operators. This fact is very important to deal with collisional plasma because the collision frequency ν must follow the Coulombian interaction law (ν ≈ 1/|v| 3 ). To perform realistic simulations in plasma physics, Coulombian interactions must be used. Therefore, relaxation operators are not relevant from a physical point of view. Moreover up to now, Asymptotic-Preserving schemes for the quasi-neutral limit have been developed either for fluid description of plasma or for collisionless plasmas.

Asymptotic-Preserving schemes have been recently used for numerous applications in the context of strong magnetic fields ( [START_REF] Degond | An Asymptotic Preserving scheme for the Euler equations in a strong magnetic field[END_REF][START_REF] Brull | Degenerate anisotropic elliptic problems and magnetised plasma simulations[END_REF]) for the gyro-fluid limit as well as in fluid mechanics for the hydrodynamic limit ( [START_REF] Jin | Efficient Asymptotic-Preserving (AP) schemes for some multi-scale kinetic equations[END_REF][START_REF] Lemou | A New Asymptotic Preserving Scheme Based on Micro-Macro Formulation for Linear Kinetic Equations in the Diffusion Limit[END_REF]). Other applications can be found for example in [START_REF] Filbet | A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Degond | An Asymptotic-Preserving method for highly anisotropic elliptic equations based on a micro-macro decomposition[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov Poisson equations[END_REF][START_REF] Brull | An asymptotic preserving scheme for a bi-fluid Euler-Lorentz system[END_REF].

Kinetic descriptions are accurate but can be too numerically expensive to be used for many real physical applications. An alternative way could be to consider a fluid description based on average quantities. Nevertheless, macroscopic descriptions are often not accurate enough. For example, in the context of inertial confinement fusion, the distribution functions considered can be far from equilibrium and in this case the fluid description is not adapted. Moreover kinetic effects like non local transport ( [START_REF] Brantov | Nonlocal heat wave propagation due to skin layer plasma heating by short laser pulses[END_REF][START_REF] Marocchino | Nonlocal heat wave propagation due to skin layer plasma heating by short laser pulses[END_REF]) or the development of some instabilities ( [START_REF] Drake | Parametric instabilities of electromagnetic waves in plasmas[END_REF]) can be important on long collisional time scales and are not captured by fluid simulations. At the same time, kinetic codes are usually limited to short time scales and cannot reach time scales studied by fluid simulations. It is therefore an important challenge to describe kinetic effects using reduced kinetic codes on fluid time scales. Then angular moments models represent intermediate models between the kinetic and fluid levels. They are less numerically expensive than kinetic models and more accurate than fluid models. They are constructed by using an angular moments extraction ( [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Minerbo | Maximum entropy Eddigton Factors[END_REF]) from the kinetic equations. But, there exists several moment models whose differences come from the choice of the closure. For example, the very popular P N closure ( [START_REF] Hauck | Positive P N closures[END_REF]) does not ensure the positivity of the distribution function. Hence we consider in this paper a M 1 moments model ( [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF]) based on an entropy minimisation principle. A M 1 -Asymptotic-Preserving scheme is also derived following the same method as in the kinetic case.

The paper is organised as follows. Section 2 introduces the Fokker-Planck-Landau-Maxwell system and its quasi-neutral limit. A reformulation of the Fokker-Planck-Landau-Maxwell system is presented in the case of one dimension in space and one dimension in velocity. The model is considered with electric fields and collision operators. Then, the method is generalised for full multi-dimensions problems with electromagnetic fields and collision operators. Section 3 introduces in detail the numerical construction of an Asymptotic-Preserving scheme for the reformulated system of section 2. Section 4 deals with the construction of an Asymptotic-Preserving scheme for the M 1 moments model from the kinetic one. Finally, section 5 presents two physically relevant numerical test cases for the M 1 -Asymptotic-Preserving scheme for different regimes. The first one shows the efficiency of the Asymptotic-Preserving scheme in the collisional quasi-neutral regime whereas the second one on the contrary corresponds to a regime where electromagnetic effects are predominant.

2. The Fokker-Planck-Maxwell system and its quasi-neutral limit.

We firstly present the kinetic model then we introduce the quasi-neutral limit.

2.1 The Fokker-Planck-Maxwell system.

Consider a kinetic model for a plasma constituted of electrons and one ion species considered fixed. Therefore the description is performed with a non-negative distribution function for electrons written f e (x, v, t) with x ∈ R n the space variable, v ∈ R n the velocity variable, n = 1, 2 or 3 and time t ≥ 0. The mass and the charge of the electron are respectively written m e and q e . The coefficient c represents the speed of light in vacuum. The Fokker-Planck-Maxwell system writes

∂f e ∂t + v.∇ x f e + q m e (E + v × B).∇ v f e = C ee (f e , f e ) + C ei (f e ), (2.1) ∂E ∂t -c 2 ∇ x × B = - j ε 0 , (2.2) ∇ x .E = q e ε 0 (n e -n i ), (2.3) ∂B ∂t + ∇ x × E = 0, (2.4) ∇ x .B = 0, (2.5)
where E and B represent respectively electric and magnetic fields. C ee and C ei are the electronelectron and electron-ion collision operators. The expression of C ee and C ei are given by

C ee (f e , f e ) = ν ee div v v ′ ∈R n S(v -v ′ )[∇ v f e (v)f e (v ′ ) -f e (v)∇ v f e (v ′ )]dv ′ , (2.6) C ei (f e ) = ν ei div v S(v)∇ v f e (v) , (2.7)
where

S(u) = 1 |u| 3 (|u| 2 Id -u ⊗ u) (2.8)
is the Landau tensor and Id is the unit tensor. The parameters ν ee and ν ei are positive physical constants.

The ionic and electronic densities n i and n e write

n i (x, t) = R n f i (x, v, t)dv, n e (x, t) = R n f e (x, v, t)dv. (2.9)
The electronic current j is given by j(x, t) = q e R n f e (x, v, t)vdv.

(2.10)

Properties of the collision operators.

The electron-electron collision operator satisfies mass, momentum and energy conservation properties (2.11) while the electron-ion collision operator satisfies only mass and energy conservation

R n C ee (f e , f e )   1 v v 2   dv = 0,
R n C ei (f e ) 1 v 2 dv = 0. (2.12)
They both dissipate the entropy i.e.

R n C ei (f e ) log f e dv ≤ 0, R n
C ee (f e , f e ) log f e dv ≤ 0, (2.13) which implies that the Boltzmann entropy Properties.

H(f e ) = R n (f e logf e -f e )dv (2.
(i) The equilibrium states of the electron-ion collision operator C ei i.e. C ei (f e ) = 0 are given by the set of isotropic functions f e (v) = f e (|v|).

(ii) The equilibrium states of the electron-electron collision operator C ee i.e. C ee (f e , f e ) = 0 are given by the Maxwellian distribution functions

f e (v) = n e ( m e 2πk B T ) 3/2 exp - m e (v -u) 2 2k B T , (2.15)
where k B is the Boltzmann constant, n e is the electronic density, T is the temperature and u represents the mean velocity.

(iii) The equilibrium states for both collision operators i.e. C ee (f e , f e ) + C ei (f e ) = 0 are given by the isotropic Maxwellian distribution function

f e (v) = n e ( m e 2πk B T ) 3/2 exp - m e v 2 2k B T . (2.16)

Scaling for the analysis of collisional processes.

For the analysis of collisional processes three important parameters are introduced: the mean free path λ ei which represents the average distance travelled by an electron between two collisions, the thermal velocity v th and the electron-ion collision frequency ν ei given by

v th = k B T m e , ν e,i = v th λ e,i . (2.17)
These parameters enable us to scale time, space and speed t = ν e,i t,

x = x/λ e,i , ṽ = v/v th . (2.18)

In the same way, we scale the electric field, the magnetic field and the distribution function

Ẽ = eE m e v th ν e,i , B = eB m e ν e,i , f = f e v 3 th n 0 . (2.19) n 0 is the initial electronic density.
With these dimensionless quantities the system (2.1) becomes the following system where we have omitted the tildes

                                       ∂f ∂t + v.∇ x f -∇ v .((E + v × B)f ) = 1 Z C e,e (f, f ) + C e,i (f ) , ∂E ∂t - 1 β 2 ∇ x × B = - j α 2 , ∂B ∂t + ∇ x × E = 0, ∇ x .E = 1 α 2 (1 -n), ∇ x .B = 0, (2.20)
where α = νe,i ωpe , ω pe represents the electronic plasma frequency, β = v th /c, n = n e /n 0 and Z the charge of the ions. In this work Z is taken equal to 1.

The dimensionless collision operators C ee (f, f ) and C ei (f ) write like in (2.6) for ν ee = ν ei = 1.

The electrostatic case.

In the electrostatic case with only one dimension for space (x ∈ R) and one for velocity (v ∈ R), the system (2.20) can be written in the following form

         ∂f ∂t + v∂ x f -E∂ v f = C e,e (f, f ) + C e,i (f ), ∂E ∂t = - j α 2 , (2.21)
where Maxwell-Poisson has to be satisfied at initial time.

Remark 1. Notice that the fourth equation of the system (2.20) called Maxwell-Gauss equation (or Poisson equation) is not used. Indeed, the second equation of (2.20) called Maxwell-Ampere equation and Poisson equation are equivalent if Poisson equation is verified at initial time.

The limit system (S 0 ) is obtained when the parameter α tends to 0 and corresponds to the quasineutral limit. It can be written in the form

     ∂f ∂t + v∂ x f -E∂ v f = C e,e (f, f ) + C e,i (f ), j = 0, (2.22)
with n = 1 at initial time.

When α tends to zero the Maxwell-Poisson equation degenerates into the algebric equation n = 1. This condition has to be satisfied at initial time.

When α is equal to zero we lose the possibility to obtain the electric field from the Maxwell-Ampere equation on collisional time scale. The limit is a singular limit, because the Maxwell-Ampere equation degenerates into an algebric equation.

Reformulation of the Maxwell-Ampere equation in the simplified case.

The aim of this part is to provide a reformulation of the Maxwell-Ampere equation that is equivalent and contains explicitly the quasi-neutral limit as a particular case when α = 0 for the electrostatic case with only one dimension for space and one for the velocity. Multiplying the first equation of (2.21) by v, integrating in velocity and using the definition of the dimensionless current

j = - R f vdv, (2.23) we obtain - ∂j ∂t + ∂ x ( R v 2 f dv) -E R v∂ v f dv = R C e,i vdv. (2.24)
In (2.11) we have seen that the calculus of R v C e,e (f, f )dv gives 0. It is important to notice that it is not the case for the electron-ion collision operator (2.12).

The derivation in time of the Maxwell-Ampere equation in the electrostatic case leads to

∂j ∂t = -α 2 ∂ 2 E ∂t 2 . (2.25)
By using (2.24), we get

α 2 ∂ 2 E ∂t 2 -E R v∂ v f dv = -∂ x ( R v 2 f )dv + R C e,i vdv. (2.26) As E R v∂ v f dv = -n e E, (2.27) the equation (2.26) becomes α 2 ∂ 2 E ∂t 2 + n e E = -∂ x ( R v 2 f )dv + R C e,i vdv. (2.28)
When the parameter α tends to 0, we find the limit problem

n e E = -∂ x ( R v 2 f dv) + R C e,i vdv. (2.29)
So the electrostatic field writes

E = -∂ x ( R v 2 f dv) + R C e,i vdv n e . ( 2

.30)

In this part we have shown that the Fokker-Planck-Maxwell system (2.21) is equivalent to the Fokker-Planck-Maxwell reformulated system

     ∂f ∂t + ∂ x (vf ) -∂ v (Ef ) = C e,e (f, f ) + C e,i (f ) , α 2 ∂ 2 E ∂t 2 + n e E = -∂ x ( R v 2 f )dv + R C e,i vdv, (2.31)
where Maxwell-Poisson has to be satisfied at initial time.

The limit system when α → 0 is the following one

     ∂f ∂t + ∂ x (vf ) -∂ v (Ef ) = C e,e (f, f ) + C e,i (f ) , E = -∂ x ( R v 2 f dv) + R C e,i vdv n e , (2.32)
where n = 1 and j = 0 at initial time.

The second equation of (2.21) imposes j = 0 when α = 0. This condition has to be satisfied at initial time.

2.6 Reformulation of the Maxwell-Ampere equation in the general case.

In this part, we generalise the method of the previous part to a non-homogeneous collisional plasma with magnetic field in multi-dimensions.

Multiplying this first equation of (2.20) by -v, integrating in velocity and using the definition of the dimensionless current (2.23) we get

- ∂j ∂t + div x ( R n v ⊗ vf dv) - R n v(E + v × B).∇ v f dv = R n C e,i (f )vdv. (2.33) As R n (v × B).∇ v f v dv = j × B, (2.34)
the same development as in the electrostatic case is performed.

The derivation in time of the Ampere-Maxwell equation in the general case leads to

∂j ∂t = -α 2 ∂ 2 E ∂t 2 + α 2 β 2 ∇ x × ∂B ∂t . (2.35)
Finally the following form is obtained

α 2 ∂ 2 E ∂t 2 + n e E -j × B = -div x ( R n v ⊗ vf dv) + α 2 β 2 ∇ x × ∂B ∂t + R n C e,i (f )vdv. (2.36)
When α tends to 0 in (2.36) we find the limit problem

n e E = -div x ( R n v ⊗ vf dv) + R n C e,i (f )vdv + j × B. (2.37)
So the electrostatic field writes

E = -div x ( R n v ⊗ vf dv) + R n C e,i (f )vdv + j × B n e . (2.38)
In this part we have shown that the Fokker-Planck-Maxwell system (2.20) is equivalent to the Fokker-Planck-Maxwell reformulated system

                     ∂f ∂t + v.∇ x f -(E + v × B).∇ v f = C e,e (f, f ) + C e,i (f ) , α 2 ∂ 2 E ∂t 2 + n e E -j × B = -div x ( R n v ⊗ vf dv) + α 2 β 2 ∇ x × ∂B ∂t + R n C e,i (f )vdv, ∂B ∂t + ∇ x × E = 0, (2.39)
where Maxwell-Poisson and Maxwell-Thomson have to be satisfied at initial time.

Remark 3. The fifth equation of the system (2.20) called Maxwell-Thomson equation is not used. Indeed the third and fifth equation of (2.20) called Maxwell-Faraday equation and Maxwell-Thomson equation are equivalent if Maxwell-Thomson equation is verified at initial time.

The limit system of (2.39) when α → 0 is the following one

           ∂f ∂t + v.∇ x f -(E + v × B).∇ v f = C e,e (f, f ) + C e,i (f ) , n e E -j × B = -div x ( R n v ⊗ vf dv) + R n C e,i vdv , ∂B ∂t + ∇ x × E = 0 , (2.40)
where n = 1 and j = 0 have to be satisfied at initial time.

The second equation of (2.40) is called the Generalised Ohm's law.

3. Discrete model.

Limitation of the classical numerical scheme.

A classical numerical scheme for the Maxwell-Ampere equation in the collisional regime writes

E n+1 i = E n i - j n i ∆t α 2 . (3.1)
The stability of this scheme depends directly on the parameter α. So, when α tends to 0, (3.1) can not be used to calculate the new electric field E n+1 i .

The aim of the following part is to establish a numerical scheme which contains explicitly the quasineutral case when α = 0. In this way, a new numerical scheme is developed for the reformulated Maxwell-Ampere equation.

3.2 Construction of an Asymptotic-Preserving Maxwell-Ampere numerical scheme.

In this part the construction of an Asymptotic-Preserving scheme for the Maxwell-Ampere reformulated equation is explained. In this first part the numerical scheme is constructed in the simplified case of an homogeneous collisionless plasma. The next part extends the method to the general non-homogeneous collisional case.

Simplified case of a homogeneous collisionless plasma.

Let us define the primal mesh M for the velocity variable v, decomposed into a family of rectangles

M p+ 1 2 =]v p-1 , v p [ ∀p ∈ {-p f , p f }
where v p = p∆v and p ∈ N represents the number of points which discretize the velocity domain. ∆v represents the energy discretisation step, which is fixed. Denote by D its associated dual mesh consisting of cells

D j =]v p-1 2 , v p+ 1 2 [ where v p-1 2 = (p - 1
2 )∆v. In the same way, a primal mesh N is defined for the space variable x, decomposed into a family of rectangles N i+ 1 2 = ]x i-1 , x i [ ∀i ∈ {1, l f } where x i = i∆x and i ∈ N represents the number of points which discretize the space domain. ∆x represents the space discretisation step, which is fixed. We denote by E its associated dual mesh consisting of cells

E i =]x i-1 2 , x i+ 1 2 [ where x i-1 2 = (i -1 2 )∆x. Let h i,p (resp. h i+ 1 2 ,p+ 1
2 ) be an

approximation of h(x i , v p ) (resp h(x i+ 1 2 , v p+ 1 2
)) for all distribution function h. Concerning the boundary conditions the velocity grid is chosen large enough to have f i,pf = f i,pf = 0 ∀i ∈ {1; lf }. It is considered that there are no particles with such velocities.

The numerical scheme for the kinetic equation is chosen in the following conservative form

f n+1 i,p -f n i,p ∆t - (E n+1 f n ) i,p+ 1 2 -(E n+1 f n ) i,p-1 2 ∆v = 0, (3.2)
where the computation of the fluxes is given by

(E n+1 f n ) i,p+ 1 2 = E n+1 i ( f n i,p + f n i,p+1 2 ) - |E n+1 i | 2 (f n i,p+1 -f n i,p ). (3.3)
Multiplying equation (3.2) by -v p ∆v and summing in p leads to

-p v p f n+1 i,p ∆v + p v p f n i,p ∆v ∆t + p v p [(E n+1 f n ) i,p+ 1 2 -(E n+1 f n ) i,p-1 2 ] = 0. (3.4)
Then using the discrete definition of the current

j i = - p v p f i,p ∆v, (3.5)
we get

j n+1 i -j n i ∆t + p v p [(E n+1 f n ) i,p+ 1 2 -(E n+1 f n ) i,p-1 2 ] = 0. (3.6)
A discrete integration by part gives

j n+1 i -j n i ∆t + p (v p -v p+1 )[(E n+1 f n ) i,p+ 1 2 ] = 0. (3.7)
By noticing that v pv p+1 = -∆v, (3.7) becomes

j n+1 i -j n i ∆t -∆v p [(E n+1 f n ) i,p+ 1 2 ] = 0. (3.8)
Using the expression (3.3) for the flux, (3.8) writes

j n+1 i -j n i ∆t -∆v p E n+1 i ( f n i,p + f n i,p+1 2 ) - |E n+1 i | 2 (f n i,p+1 -f n i,p ) = 0. (3.9) Remark 3. It is important to notice that p |E n+1 i | 2 (f n i,p+1 -f n i,p ) = |E n+1 i | 2 p (f n i,p+1 -f n i,p ) = 0. (3.10)
Therefore, no linearisation nor approximation is required to compute E n+1 i from (3.9). Using (3.10), (3.9) becomes

j n+1 i -j n i ∆t -∆vE n+1 i p ( f n i,p + f n i,p+1 2 ) = 0. (3.11)
The following scheme for the Maxwell-Ampere equation is used

E n+1 i -E n i ∆t = - j n+1 i α 2 . (3.12)
Contrarily to the classical scheme (3.1) the current j in (3.12) is chosen implicit.

So from (3.11), it comes that

-α 2 E n+1 i -2E n i + E n-1 i ∆t 2 -∆vE n+1 i p ( f n i,p + f n i,p+1
2 ) = 0, (3.13) which is the numerical scheme for the reformulated Maxwell-Ampere equation in the simplified case of a homogeneous collisionless plasma.

By letting α to 0 in (3.13) we recover

E n+1 i = 0. (3.14)

Generalisation to a non-homogeneous collisional plasma.

In this part the method is generalised to the case of a non-homogeneous collisional plasma. An Asymptotic-Preserving scheme is constructed for the second equation of (2.31).

By using a conservative discretisation for the Fokker-Planck equation we obtain

f n+1 i,p -f n i,p ∆t + (vf n ) i+ 1 2 ,p -(vf n ) i-1 2 ,p ∆x - (E n+1 f n ) i,p+ 1 2 -(E n+1 f n ) i,p-1 2 ∆v = C n ee,i + C n ei,i , (3.15)
where

C n ei,i = 1 ∆v S p+ 1 2 f n i,p+1 -f n i,p ∆v -S p-1 2 f n i,p -f n i,p-1 ∆v (3.16)
and

S p+ 1 2 = S( v p + v p+1 ∆v ). (3.17)
The expression of S is given in (2.8). The numerical scheme for the operator C ee,i is not given, because using (2.11) this term cancels in the calculation.

Using upwind fluxes, (3.15) leads to

f n+1 i,p -f n i,p ∆t + v p f n i+1,p -f n i-1,p -|v p | f n i+1,p -2f n i,p + f n i-1,p 2∆x - E n+1 i f n i,p+1 -f n i,p-1 -|E n+1 i | f n i,p+1 -2f n i,p + f n i,p-1 2∆v = C n ee,i + C n ei,i .
Multiplying the previous equation by -v p ∆v and summing in p leads to

-p v p f n+1 i,p ∆v + p v p f n i,p ∆v ∆t - ∆v 2∆x p v 2 p f n i+1,p -f n i-1,p -v p |v p | f n i+1,p -2f n i,p + f n i-1,p + 1 2 p v p E n+1 i f n i,p+1 -f n i,p-1 -|E n+1 i |v p f n i,p+1 -2f n i,p + f n i,p-1 (3.18) = - p C n ei,i v p ∆v.
Using (3.5), the computation of the previous equation leads to

j n+1 i -j n i ∆t - ∆v 2∆x p v 2 p f n i+1,p -f n i-1,p -v p |v p | f n i+1,p -2f n i,p + f n i-1,p + 1 2 p v p E n+1 i f n i,p+1 -f n i,p-1 -|E n+1 i |v p f n i,p+1 -2f n 1,i,p + f n i,p-1 (3.19) = - p C n ei,i v p ∆v.
By using (3.12), we get

-α 2 E n+1 i -2E n i + E n-1 i ∆t 2 - ∆v 2∆x p v 2 p f n i+1,p -f n i-1,p -v p |v p | f n i+1,p -2f n i,p + f n i-1,p + 1 2 p v p E n+1 i f n i,p+1 -f n i,p-1 -|E n+1 i |v p f n i,p+1 -2f n 1,i,p + f n i,p-1 (3.20) = - p C n ei,i v p ∆v. Moreover p |E n+1 i |v p f n i,p+1 -2f n i,p + f n i,p-1 = 0. (3.21)
Indeed a discrete integration by part gives

p |E n+1 i |v p f n i,p+1 -2f n i,p + f n i,p-1 = |E n+1 i | p v p (f n i,p+1 -f n i,p ) - p v p+1 (f n i,p+1 -f n i,p ) , = |E n+1 i | p (v p -v p+1 )(f n i,p+1 -f n i,p ) , = -|E n+1 i |∆v p (f n i,p+1 -f n i,p ) , = 0.
Finally, the Asymptotic-Preserving scheme for the second equation of (2.31) writes

-α 2 E n+1 i -2E n i + E n-1 i ∆t 2 - ∆v 2∆x p v 2 p f n i+1,p -f n i-1,p -v p |v p | f n i+1,p -2f n i,p + f n i-1,p + E n+1 i 2 p v p f n i,p+1 -f n i,p-1 = - p C n ei,i v p ∆v.
In the limit case when α tends to zero, the scheme becomes

E n+1 i = ∆v ∆x p v 2 p f n i+1,p -f n i-1,p -v p |v p | f n i+1,p -2f n i,p + f n i-1,p -2 p C n ei,i v p ∆v p v p f n i,p+1 -f n i,p-1 . (3.22)
In the case of an homogeneous collisionless plasma the expression (3.14) is recovered.

Generalisation to an non-homogeneous collisional plasma with electromagnetic fields.

In this part we derive the numerical scheme for the reformulated Maxwell-Ampere equation in the simplified case of 1 dimension in space and 3 dimensions in velocity. The scheme can be extended to the case of 3 dimensions in space. We consider a cartesian case with an electric and a magnetic field of the form E = (E x (t, x, y), E y (t, x, y), 0), B = (0, 0, B z (t, x, y)). (3.23) Following the same method as for the electrostatic case, we derive the following numerical scheme for the reformulated Maxwell-Ampere equation

-α 2 E n+1 x,l -2E n x,l + E n-1 x,l ∆t 2 -∆v x ∆v y ∆v z i,j,k E n+1 x,l + j∆v y B n+1 z,l f n l,i,j,p -∆v 2 x ∆v z i,j,k v x,i E n+1 y,l -i∆v x B n+1 z,l f n l,i,j,p = - i,j,k C n ei,i v x,i ∆v x ∆v y ∆v z + ∆v x ∆v y ∆v z 2∆x i,j,k v 2 x,i f n l+1,i,j,k -f n l-1,i,j,k -v x,i |v x,i | f n l+1,j,k,p -2f n l,i,j,k + f n l-1,i,j,k , -α 2 E n+1 y,l -2E n y,l + E n-1 y,l ∆t 2 -∆v 2 y ∆v z i,j,k E n+1 x,l + j∆v y B n+1 z,l f n l,i,j,p -∆v x ∆v y ∆v z i,j,k v x,i E n+1 y,l -i∆v x B n+1 z,l f n l,i,j,p = - i,j,k C n ei,i v y,i ∆v x ∆v y ∆v z + α 2 β 2 ∆t B n+1 z,l+1 -B n+1 z,l-1 2∆x - B n z,l+1 -B n z,l-1 2∆x + ∆v x ∆v y ∆v z 2∆x i,j,k v y,i v x,i f n l+1,i,j,k -f n l-1,i,j,k -v y,i |v x,i | f n l+1,i,j,k -2f n l,i,j,k + f n l-1,i,j,k
, where l is the index for space, i the index for the first coordinate in speed, j for the second and k for the third. Also ∆t, ∆x, ∆v x , ∆v y , ∆v z are respectively the time step, the space step, the velocity step in the first, second and third dimension. In this case there are two equations, we notice they are coupled.

4. The Asymptotic-Preserving scheme for the M 1 -Maxwell moment model. This part is devoted to the derivation of an Asymptotic-Preserving scheme for the M 1 model associated to the system (2.20). For the sake of clarity, we firstly recall the derivation of the M 1 model that is performed in ( [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Fokker-Planck-Landeau equation of electrons[END_REF][START_REF] Mallet | General moment system for plasma physics based on minimum entropy principle[END_REF]) and next we present the numerical scheme.

f 0 (ζ) = ζ 2 1 -1 f (µ, ζ)dµ, f 1 (ζ) = ζ 2 1 -1 f (µ, ζ)µdµ, f 2 (ζ) = ζ 2 1 -1 f (µ, ζ)µ 2 dµ. (4.1)
For moments model one fundamental point is the definition of the closure which guarantees that the highest moment writes as a function of the previous ones. This closure corresponds to an approximation of the distribution function from which the moments system is constructed. In the present paper the M 1 model closure is obtained from an entropy minimisation principle.

Closure for the M 1 moment model.

In order to close the problem we need to define f 2 as a function of f 0 and f 1 . For the M 1 model, the construction of the closure is based on an entropy minimum principle ( [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Minerbo | Maximum entropy Eddigton Factors[END_REF]). Indeed the closure is obtained by solving

min g≥0 { H(g) / ∀ζ ∈ R + , 1 -1 g(µ, ζ)dµ = f 0 (ζ), 1 -1 g(µ, ζ)µdµ = f 1 (ζ) }, (4.2)
where H(g) is the Boltzmann entropy defined in (2.14).

The entropy minimum principle implies that the solution of (4.2) writes ( [START_REF] Dubroca | Étude théorique et numérique d'une hiéarchie de modèles aux moments pour le transfert radiatif[END_REF][START_REF] Mallet | General moment system for plasma physics based on minimum entropy principle[END_REF])

f (µ, ζ) = ρ(ζ)exp(-µ a 1 (ζ)), (4.3)
where ρ(ζ) is a positive scalar, and a 1 (ζ) a real valued scalar.

Referring to [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF], f 2 writes

f 2 = χ(f 0 , f 1 )f 0 with χ = |a 1 | 2 -2|a 1 |coth(|a 1 |) + 2 |a 1 | 2 .

Collision operators.

Moments extraction of the electron-electron collision operator (2.6) is complicated due to its nonlinearity. That is why, it is often simplified by considering only the isotropic part of the distribution function ( [START_REF] Buet | Conservative and entropy decaying numerical scheme for the isotropic Fokker-Planck-Landau equation[END_REF]). Nevertheless, it has been shown ( [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Fokker-Planck-Landeau equation of electrons[END_REF]) that this simplification does not preserve the realisibility domain which guarantees that the solution of the moments system can be expressed as the moments of a nonnegative distribution function. Consequently, the collisional operators that are used ( [START_REF] Mallet | General moment system for plasma physics based on minimum entropy principle[END_REF][START_REF] Mallet | An entropic scheme for an angular moment model for the classical Fokker-Planck-Landeau equation of electrons[END_REF]), are constructed according to a linearisation around the equilibrium state of C ei . So C ee is approached by Q ee which correspond to its linearisation around the isotropic part of the distribution function Q ee .

Hence, Q ee writes

Q e,e (f ) = 1 ζ 2 ∂ ζ ζ +∞ 0 J(ζ, ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ f (ζ) -f (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ζ ′2 dζ ′ , (4.4) with J(ζ, ζ ′ ) = 2 3 inf ( 1 ζ 3 , 1 ζ ′3 )ζ ′2 ζ 2 . (4.5)
For the electron-ion operator, no approximation is performed because this operator is already linear. In the following Q e,i will replace C e,i to be consistent with the notation (4.4) of Q e,e .

Q ei (f ) = 1 ζ 3 ∂ ∂µ (1 -µ 2 ) ∂f ∂µ . (4.6)
4. [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetised plasma simulations[END_REF] The M 1 moment model.

The angular integration ([27]) leads to

   ∂ t f 0 + ∇ x .(ζf 1 ) -∂ ζ (Ef 1 ) = Q 0 (f 0 ), ∂ t f 1 + ∇ x .(ζf 2 ) -∂ ζ (Ef 2 ) + E (f 0 -f 2 ) ζ = Q 1 (f 1 ) + Q 0 (f 1 ), (4.7)
where the collisional operators Q 0 and Q 1 are given by

Q 0 (f 0 ) = 2 3 ∂ ζ ζ 2 A(ζ)∂ ζ ( f 0 ζ 2 ) -ζB(ζ)f 0 , (4.8) Q 1 (f 1 ) = - 2f 1 ζ 3 . (4.9)
The coefficients A(ζ) and B(ζ) write In this part the reformulation of the Maxwell-Ampere equation for the M 1 model is detailed. Considering a conservative scheme for the second equation of (4.7) we write

A(ζ) = ∞ 0 min( 1 ζ 3 , 1 ω 3 )ω 2 f 0 (ω)dω, (4.10) B(ζ) = ∞ 0 min( 1 ζ 3 , 1 ω 3 )ω 3 ∂ ω ( f 0 (ω) ω 2 )dω. ( 4
f n+1 1,i,p -f n 1,i,p ∆t + (ζf n 2 ) i+ 1 2 ,p -(ζf n 2 ) i-1 2 ,p ∆x - (E n+1 f n+1 2 ) i,p+ 1 2 -(E n+1 f n+1 2 ) i,p-1 2 ∆ζ (4.12) + E n+1 ζ p (f n+1 0,i,p -f n+1 2,i,p ) = Q n 1,i,p + Q n 0,i,p .
The discrete collision operators involved in (4.12) are respectively given by

Q n 1,i,p =- 2f n 1,i,p ζ 3 p , Q n 0,i,p = 2 3∆ζ p (ζ 2 p+ 1 2 A(ζ p+ 1 2 ) 1 ∆ζ p+ 1 2 f n 1,i,p+1 ζ 2 p+1 - f n 1,i,p ζ 2 p -ζ p+ 1 2 B(ζ p+ 1 2 )f n 1,i,p+ 1 2
)

-(ζ 2 p-1 2 A(ζ p-1 2 ) 1 ∆ζ p-1 2 f n 1,i,p ζ 2 p - f n 1,i,p-1 ζ 2 p-1 -ζ p-1 2 B(ζ p-1 2 )f n 1,i,p- 1 2
) .

Using a HLL scheme for the fluxes in (4.12), it holds that

f n+1 1,i,p -f n 1,i,p ∆t + ζ p f n 2,i+1,p -f n 2,i-1,p -|ζ p | f n 1,i+1,p -2f n 1,i,p + f n 1,i-1,p 2∆x - E n+1 i f n 2,i,p+1 -f n 2,i,p-1 -|E n+1 i | f n 1,i,p+1 -2f n 1,i,p + f n 1,i,p-1 2∆ζ (4.13) + E n+1 i ζ p (f n 0,i,p -f n 2,i,p ) = Q n 1,i,p + Q n 0,i,p .
Multiplying the previous equation by -ζ p ∆ζ and summing in p leads to

-p ζ p f n+1 1,i ∆ζ + p ζ p f n 1,i ∆ζ ∆t - ∆ζ 2∆x p ζ 2 p f n 2,i+1,p -f n 2,i-1,p -ζ 2 p f n 1,i+1,p -2f n 1,i,p + f n 1,i-1,p + (4.14) 1 2 p ζ p E n+1 i f n 2,i,p+1 -f n 2,i,p-1 -|E n+1 i |ζ p f n 1,i,p+1 -2f n 1,i,p + f n 1,i,p-1 - p E n+1 i (f n 0,i,p -f n 2,i,p )∆ζ = - p ζ p Q n 1,i,p ∆ζ.
Here again, the term containing the electron-electron collision operator cancels.

We use the definition of the dimensionless current

j j = - R + f 1 ζdζ, (4.15)
which can be written on the discrete form

j n i = - p f n 1,i,p ζ p ∆ζ. (4.16)
Therefore the scheme (4.14) becomes

j n+1 i -j n i ∆t - ∆ζ 2∆x p ζ 2 p f n 2,i+1,p -f n 2,i-1,p -ζ 2 p f n 1,i+1,p -2f n 1,i,p + f n 1,i-1,p + 1 2 p ζ p E n+1 i f n 2,i,p+1 -f n 2,i,p-1 -|E n+1 i |ζ p f n 1,i,p+1 -2f n 1,i,p + f n 1,i,p-1 (4.17) - p E n+1 i (f n 0,i,p -f n 2,i,p )∆ζ = - p ζ p Q n 1,i,p ∆ζ.
Using the scheme (3.12), the expression (4.17) becomes

-α 2 E n+1 i -2E n i + E n-1 i ∆t 2 - ∆ζ 2∆x p ζ 2 p f n 2,i+1,p -f n 2,i-1,p -ζ 2 p f n 1,i+1,p -2f n 1,i,p + f n 1,i-1,p + 1 2 p ζ p E n+1 i f n 2,i,p+1 -f n 2,i,p-1 -|E n+1 i |ζ p f n 1,i,p+1 -2f n 1,i,p + f n 1,i,p-1 (4.18) -E n+1 i p (f n 0,i,p -f n 2,i,p )∆ζ = - p ζ p Q n 1,i,p ∆ζ.
Like for the kinetic scheme in (3.21), it holds that

p |E n+1 i |ζ p f n 1,i,p+1 -2f n 1,i,p + f n 1,i,p-1 = 0. (4.19)
Therefore the final scheme obtained reads

E n+1 i = -α 2 (2E n i -E n-1 i ) ∆t 2 + β 1 (f n 0,i , f n 1,i ) - α 2 ∆t 2 + β 2 (f n 0,i , f n 1,i ) , (4.20)
where the coefficients β 1 and β 2 are given by

β 1 = ∆ζ 2∆x p ζ 2 p f n 2,i+1,p -f n 2,i-1,p -ζ 2 p f n 1,i+1,p -2f n 1,i,p + f n 1,i-1,p - p ζ p Q n 1,i,p ∆ζ, (4.21) β 2 = 1 2 p ζ p f n 2,i,p+1 -f n 2,i,p-1 - p (f n 0,p,i -f n 2,p,i )∆ζ. (4.22)
Remark 5. The stability of this new scheme does not depend on the parameter α. So, the electrostatic field can be obtained even if α becomes equal to zero. Remark 6. Following the same procedure as for the Fokker-Planck-Maxwell system, this reformulation can be generalised for multi-dimension problems with magnetic fields.

Numerical test cases.

This section presents two physically relevant numerical experiments where opposite regimes are considered. The first one studies the relaxation of a localised temperature profile in the quasi-neutral regime, in this regime collisions between particles are predominant. The second one studies the two beams interaction, in this case electrostatic effects are predominant. These numerical experiments show the efficiency of the new scheme for different regimes.

Hot spot relaxation.

We study the relaxation of a localised temperature profile. This phenomenon investigated by Batishchev ([1], [START_REF] Brantov | Nonlocal heat wave propagation due to skin layer plasma heating by short laser pulses[END_REF]) corresponds physically to the heating of a plasma during a short time and to the relaxation phenomenon which follows. The important temperature gradients due to the localised heating induce a non-local heat transport. Here, we consider the collisional regime. This configuration is particularly interesting because it enables to study the coupling of the M 1 model with the Maxwell-Ampere Asymptotic-Preserving scheme.

Initially the distribution function for electrons is a Maxwellian with a Gaussian temperature profile

T e (x, t = 0) = T 0 + T 1 exp(- x 2 D 2 ), (5.1)
with T 0 = 1Kev, T 1 = 4Kev and D = 8, 44λ ei . There is no electric field at the initial time. Here we choose specular reflection conditions for boundary conditions. The space scale goes from -80λ ei to 80λ ei . The energy space scale goes from 0 to 50 v th .

The following plots show the evolution of the temperature and electric field profiles until 30 τ ei . Then at t = 2 τ ei , we observe that the temperature profile starts to relax to a colder temperature and the electric field which is proportional to the gradient of temperature also decreases. The numerical scheme represents the good behaviour of the hot spot relaxation phenomenon. The results of our M 1 -Asymptotic-Preserving scheme (M 1 -AP) have been compared with the ones obtained by a kinetic code ( [START_REF] Marocchino | Nonlocal heat wave propagation due to skin layer plasma heating by short laser pulses[END_REF]). On the following plots, the temperature and the electrostatic field profiles are represented as function of space for different times. The (M 1 -AP) results are given in green while the kinetic results are in red. Both results show a good agreement. Small differences are observed concerning the amplitude of the temperature and the electric field. The relaxation phenomenon observed with the (M 1 -AP) code is faster than the one with the kinetic code. It is interesting to notice that there is a large calculation time difference. The simulation with the kinetic code requires the use of 50 processes during several days while the (M 1 -AP) code only needs few minutes with one process.

Thanks to the rapidity of the M 1 Asymptotic-Preserving code a mesh convergence study has been performed. With 500 points for the space grid and 80 points for the energy grid the results are converged. The time step used is ∆t = 10 -3 τ ei in order to respect classical stability conditions. Asymptotic-Preserving scheme.

Remark 7. In this case the parameter α which represents the ratio between the electron-ion collision frequency and electronic plasma frequency is equal to 4.10 -4 . In order to avoid a severe constraint on the time step we use the new M1-Asymptotic-Preserving scheme. With the same CFL conditions, the classic Maxwell-Ampere numerical scheme breaks down from the very first iterations. Remark 8. It is important to notice that the Asymptotic-Preserving scheme is stable even in the case α = 0.

Two electron beams interaction.

In this part we study the interaction between two electron beams. This collisionless test case enables us to study the regime where electrostatic effects are predominant. Therefore for this test case we have C ee = C ei = 0. Consider two electron beams propagating at velocity v 0 and v 1 . The dispersion relation is given by 1

- 1 (ω -kv 0 ) 2 - 1 (ω -kv 1 ) 2 = 0, (5.2)
where v 0 and v 1 denote the beams velocities. This configuration can lead to electrostatic instabilities. Indeed, the solutions of the form Ae iωt+ikx are unstable when ω I the imaginary part of ω is strictly positive. In the case v 0 = -v 1 we can show that the solution is stable if

kv 0 ≥ √ 2.
This test is problematic for the M 1 model. Indeed, if we consider two electron beams propagating with opposite velocities the distribution function is well defined. Nevertheless, the M 1 model considers only the angular moments f 0 and f 1 . For the calculus of f 1 the two populations contributions cancel and we get f 1 = 0. The M 1 model sees an isotropic configuration which is not the case. To overcome this problem we use the superposition principle that is valid because the model is linear. Two particle populations (one per beam) are considered. For each time step the M 1 problem is solved for the first population then for the second one. Hence the Maxwell equations are solved taking into account the two distribution functions.

In the case of two streams propagating with opposite velocities v d and -v d , the initial conditions are f (x, v, t = 0) = 0. The parameter A is introduced to perturb the initial condition in order to enable the development of the electrostatic instability. The energy scales goes from 0 to 12 v th and the space scale from 0 to 25 ω pe . With 100 points for the space grid and 128 points for the energy grid the results are converged. On the following plots the distribution function is represented in the phase space for the initial time and the final time t = 30 plasma periods. In this example v d = 4, A = 0.001 and periodical boundary conditions are used. On the second plot the interaction between the two streams is observed. Our results have been compared with a kinetic code ( [START_REF] Duclous | High order resolution of the Maxwell-Fokker-Planck-Landau model intended for ICF applications[END_REF]). On Figure 4, the evolution of the electrostatic energy is represented as a function of time for the (M 1 -AP) code in green and for a kinetic code in red. The first plot shows the results for A=0.001 and the second one for A=0.1.

In the case of small perturbations (A=0.001), the M 1 model and the kinetic code give analogous results. In the case of strong perturbations (A=0.1), the (M 1 -AP) code and the kinetic code show some differences after a long time. In the case of a strong perturbation, a non-linear regime is obtained and it is well-known that the M 1 model is not accurate enough ( [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF]). This numerical experiment shows the good behaviour of the (M 1 -AP) scheme in a regime where electrostatic effects are predominant.

Conclusion.

In this work, we have constructed an Asymptotic-Preserving scheme for the full Fokker-Planck-Landau-Maxwell system which handles the quasi-neutral limit without any contraction of time and space steps. Remark that this model is considered with the real collisional operators. This fact is important in plasma physics because the model is relevant for Coulombian interactions. We have first established a reformulated Fokker-Planck-Landau-Maxwell system then used it to construct the Asymptotic-Preserving scheme. The method has been extended to the general case of collisional plasmas in electromagnetic fields for multi-dimensions problems. An M 1 -Asymptotic-Preserving scheme has been derived. Next, the (M 1 -AP) scheme has been implemented and two numerical test cases have been performed. The first one shows the efficiency of the Asymptotic-Preserving scheme in the quasi-neutral regime. The second one on the contrary corresponds to a regime where electromagnetic effects are predominant. The scheme, accurate and fast, works in both regimes.

In this paper and in previous ones ( [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Fokker-Planck-Landeau equation of electrons[END_REF][START_REF] Mallet | General moment system for plasma physics based on minimum entropy principle[END_REF]), the ions have been considered as fixed due to their important mass compared to the ones of electrons. The model must be extended when times intervals during which their motion can not be neglected are considered. Therefore, we expect to develop in a forthcoming paper a model taking into account ions motion.

  [START_REF] Degond | Numerical approximation of the Euler-Maxwell model in the quasineutral limit[END_REF] is a Lyapunov function for the equation (2.1).

4. 1

 1 Moment models. If S 2 is the unit sphere, Ω = v |v| represents the direction of propagation of the particle and µ = Ω x = cos θ, θ ∈ [0, π]. In this part, we choose a one dimensional direction of propagation, i.e we take µ ∈ [-1, 1] as the direction of propagation instead of Ω. By setting, ζ = |v| the distribution function f writes in spherical coordinates f (µ, ζ, x). Hence the three first angular moments of the distribution function are given by

  .11) 4.5 A numerical scheme for the M 1 model.

Figure 1 .

 1 Figure 1. Representation of the temperature and electric field as function of space for different times.

Figure 2 .

 2 Figure 2. Comparison of the temperature and electric field for a kinetic code and the M 1Asymptotic-Preserving scheme.

5 [ 1 +

 51 Acos(kx)M v d (v) + (1 -Acos(kx))M -v d (v)], (5.3) with M ±v d (v) = n e ( m e 2πk B T e ) 3/2 exp -m e (v ∓ v d ) 2 2k B T e . (5.4)

Figure 3 .

 3 Figure 3. Distribution function as function of space and velocity at initial time (left) and after 30 plasma periods (right).

Figure 4 .

 4 Figure 4. Temporal evolution of the electrostatic energy in the linear regime (top) and in the non-linear regime (bottom).

Asymptotic-Preserving scheme for the Fokker-Planck-Landau-Maxwell system in the quasi-neutral regime.

Acknowledgements. We are thankful with Vladimir Tikhonchuk, Axel Klar and Cory Hauck for fruitful discussions.