Sylvain Crovisier
email: sylvain.crovisier@u-bourgogne.fr

Micha L Rams
email: rams@impan.gov.pl

IFS attractors and Cantor sets

Keywords:

We build a metric space which is homeomorphic to a Cantor set but cannot be realized as the attractor of an iterated function system. We give also an example of a Cantor set K in R 3 such that every homeomorphism f of R 3 which preserves K coincides with the identity on K.

Introduction

We are interested in the problem of characterization of compact sets that are limit sets of hyperbolic dynamical systems. In this paper we study iterated function systems that are simplified models for the smooth hyperbolic dynamics. Previous works (see [DH, K]) investigated which compact metric spaces can be attractors of iterated function systems on some euclidean space. We would like to carry on this discussion with the following question:

Is every Cantor set an attractor of some iterated function system? Let us first recall that a continuous mapping f of a metric space (X, d) into itself is a contractive map if there exists a constant σ ∈ (0, 1) such that for any points x, y in X, the distance d(f (x), f (y)) is less or equal to σd(x, y). An iterated function system (or IFS) is a finite family of contractive maps {f 1 , . . . , f s } acting on a complete metric space. It is well known (see [H]) that such a dynamics possesses a unique compact set K ⊂ X which is non-empty and invariant by the IFS, in the following sense:

K = ∪ s i=1 f i (K).
1 One calls this compact set the limit set or the attractor of the IFS.

As an example, the middle one third Cantor set is the attractor of the IFS {f 1 , f 2 } on R defined by f 1 (x) = x/3, f 2 (x) = (x + 2)/3.

Our first result shows that for some other metrics, the Cantor set is no more the attractor of an IFS.

Theorem 1.1. There exists a Cantor set X 1 and a Borel probability measure µ supported on X 1 such that for every contractive map f : X 1 → X 1 , we have µ(f (X 1 )) = 0. This set X 1 cannot be an attractor of an iterated function system {f 1 , . . . , f s } (even if one allows countably many maps in the definition of the IFS) since X 1 has full measure by µ but f 1 (X 1 ) ∪ • • • ∪ f s (X 1 ) has zero measure by µ. The example may be generalized:

Corollary 1.2. For any iterated function system on a complete metric space (X, d), the attractor K is not isometric to the Cantor set X 1 , defined in Theorem 1.1.

In the previous case, the obstruction was metrical. If one specifies the Cantor set K and the ambient space X there may exist also topological obstructions for K to be an attractor of an IFS defined on X, even if X is a smooth manifold such as R d .

Theorem 1.3. There exists a Cantor set

X 2 ⊂ R 3 such that if f is a home- omorphism of R 3 which satisfies f (X 2 ) ⊂ X 2 then f |X 2 = id.
In particular a finite set of homeomorphisms of R 3 can not be an IFS whose attractor is X 2 . The set X 2 here is a variation on Antoine's necklace.

This example can be easily generalized to higher dimensions but in dimension one or two the situation is completely different: Proposition 1.4. For any Cantor set X ∈ R (or in R 2 ) and any two points x, y ∈ X, there exists a contractive homeomorphism f : R → R (or

f : R 2 → R 2 ) such that f (X) ⊂ X and f (x) = y.

Constructions

Given Y , subset of a metric space X we denote its complement by (Y ) c , its interior Int(Y ), its boundary ∂(Y ) and (in the case Y is bounded) its diameter diam(Y ). We will denote by B(z, r) the open ball centered at z ∈ X with radius r.

g 1 g 2 X1 I (2) I (1) I (0) Figure 1: Construction of X 1 2.1 Proof of Theorem 1.1 2.
1.1 Definition of the Cantor set X 1 and the measure µ

The Cantor set X 1 is obtained as the intersection of a decreasing sequence (I (k) ) of compact sets in R. Each set I (k) is a finite union of pairwise disjoint compact intervals I (k) i that have the same length. We will construct inductively the sequence (I (k) ). The first set

I (0) = [0, 1]. Given I (k) = I (k)
j ,we choose n k+1 closed intervals I (k+1) i inside every I (k) j . We demand those intervals to be pairwise disjoint and of equal length, that their union contains the endpoints of I (k) j and also that the gaps between them are of equal length g k+1 . The set I (k+1) will be the union of all I (k+1) i . Obviously, the bounded components of R \ I (k+1) are either gaps created at the previous steps of the construction or new gaps of size g k+1 each.

We take care along the induction that in each consecutive step the number of intervals n k increases while the size of gaps g k decreases. We define X 1 = k I (k) which is obviously a Cantor subset of R, see the figure 1. We define the measure µ as to be equidistributed: all the intervals I (k) i of a given level k have the same measure, equal to the reciprocal of the total number of intervals of that level, ie. µ(I

(k) i ) = k j=1 n -1 j .
By the Kolmogorov theorem ( [B2]) it uniquely defines a probability measure supported on X 1 . 3 2.1.2 Measure of f (X 1 ) Now, we prove that X 1 satisfies the assertion of Theorem 1.1: let f be a contraction from X 1 into itself and let

X (k) i = X 1 ∩ I (k) i be the part of X 1 contained in one of the k-th level intervals I (k) i . We claim that f (X (k) i ) must be contained in some (k + 1)-th level interval.
Assuming it is not the case, f (X

(k)
i ) intersects at least two (k + 1)-th level intervals. Since (g n ) is decreasing, these intervals are in distance at least g k+1 from each other. Hence, f (X (k) i ) may be divided onto two subsets A, B such that the distance between any point from A and any point from B is not smaller than g k+1 . As the map f is contracting, the preimages of A and B (covering together whole X (k) i ) must lie in distance strictly greater than g k+1 . But if we could divide X (k) i into two such subsets, this would imply the existence of a gap inside X (k) i of size strictly greater than g k+1 . By construction, such a gap would be created at a step ℓ larger or equal to k + 1 and would be of length g ℓ . The sequence (g n ) would not be strictly decreasing. This contradiction proves the claim.

Hence, as there are k j=1 n j intervals of level k in I (k) , the whole set f (X 1 ) is contained in the union of at most k j=1 n j (k + 1)-th level intervals. This implies

µ(f (X 1 )) ≤ k i=1 n i k+1 j=1 n -1 j = n -1 k+1
and the right hand side tends to zero when k tends to infinity.

Proof of Theorem

1.3 2.2.1 Definition of the Cantor set X 2
The construction is a variation on Antoine's necklace example built in [A 1 , A 2 ] (see also for example [M], chapter 18 or [B1], section IV.7 for more details).

We start from a filled compact torus T (0) ⊂ R 3 that will also be noted T ∅ . In the first step we find in the interior of T (0) some number (say, n = n ∅ > 2) of disjoint compact tori T 1 , . . . , T n linked together to form a closed chain going around the torus T (0) , see figure 2. They will be called first level tori. We denote their union by T (1) . In the second step we take inside each of the first level tori T i another chain of smaller tori (called second level tori) going around (of n 1 elements inside T 1 , n 2 inside T 2 and so on up to n n ∅ ). They will be denoted by T i,1 , . . . , T i,n i and their union by T (2) ⊂ T (1) .

Figure 2: The first sets T (0) and T (1) (here we have n ∅ = 4)

We repeat the procedure inductively, taking care that the diameters of the tori we use in the construction go down to 0 and that, at any step, the lengths of the chains are greater than 2, different from each other, and different from the lengths of all other chains introduced at the previous steps (in Antoine's example one takes chains of length four at every step). All the tori of level k are described by finite sequences ω k ∈ Σ of length k where:

Σ = k∈N {ω k = (ω 1 , . . . , ω k ) ∈ N k ; ∀ j≤k 1 ≤ ω j ≤ n ω k-1 }.
(We denoted by ω k a sequence of length k, by ω j the j-th element of the sequence and by ω j the subsequence formed by the j first terms in ω k .)

We define

X 2 = k T (k) .
As one can easily check from the construction, X 2 is a Cantor set. We introduce also the rings of level k in X 2 defined by

Y ω k = X 2 ∩ T ω k .
Let us remark that the sets Y ω k are open and closed in X 2 . They define a basis for the topology on X 2 . A chain C in X 2 is a union of rings of same level:

C = Y ω k ,i 1 ∪ Y ω k ,i 2 ∪ • • • ∪ Y ω k ,ir
where ω k belongs to Σ and where {i 1 , . . . , i r } is an interval in Z/n ω k Z.

Topological properties of X 2

Let Y be a compact set of R 3 . We will say that Y may be cleaved if there exists a decomposition Y = A 1 ∪ A 2 of Y in two compact disjoint and non empty sets and some isotopy of homeomorphisms (h t ) t∈[0,1] of such that h 0 = id and such that h 1 (A 1 ) and h 1 (A 2 ) are contained in two disjoint euclidean balls of R 3 . We prove in this section the following topological characterization of rings of X 2 :

Proposition 2.1. A compact subset Y of X 2 is a ring of X 2 if and only if it satisfies the following properties: a) Y can not be cleaved. b) Y is "cyclic": there exists a partition of Y in three disjoint compact non empty subsets Y = A 1 ∪ A 2 ∪ A 3 such that each A i and each union A i ∪ A j can not be cleaved.
We first check that rings satisfy these properties.

Lemma 2.2. No chain C of X 2 may be cleaved.

In particular the rings can not be cleaved.

Proof. In [A 1 , A 2 ] (see also [M], problem 18.2), Antoine proved that any 2sphere in R 3 that do not intersect the Antoine's necklace can not separate two points of the necklace from one another. By the same argument, this property is also satisfied by any chain of X 2 . This implies the lemma.

We then check that rings are "cyclic" (i.e. satisfy Proposition 2.1 b)): let Y ω k be a ring of X 2 . We set

A 1 = Y ω k ,1 , A 2 = Y ω k ,2 and A 3 = Y ω k \ (A 1 ∪ A 2 ).
All the sets A i or A i ∪ A j are chains of X 2 and can not be not cleaved from lemma 2.2, what was to be shown.

The other part of the proof of Proposition 2.1 uses the following lemma.

Lemma 2.3. Let Y be a compact and proper subset of some ring Y ω k in X 2 . Then, there exists an isotopy (h t ) t∈[0,1] in the space of homeomorphisms such that h t coincides with the identity on (Int(T ω )) c , h 0 = id and h 1 sends Y on a small euclidean ball included in Int(T ω ).

Proof. This lemma is clear when Y is a chain Y ω k ,i 1 , . . . , Y ω k ,ir of T ω k . In the general case, we note that Y ω k \ Y is open in Y ω k and contained in some ring. Hence, Y is contained in a finite union of pairwise disjoint rings that is strictly included in Y ω k . One observes that it is enough to prove the assertion for such union of rings, it will imply the result for Y . Note that such a union of rings is a finite union of pairwise disjoint chains not linked to each other. The proof is then done by induction on the number of chains : from the remark we did at the beginning of this proof, one starts by shrinking the chains of largest order inside small disjoint balls. This allows us to consider the chains of next largest order and to shrink them. Repeating this procedure, one concludes the proof.

We get a counterpart to lemma 2.2 Lemma 2.4. Let Y be a closed subset of X 2 . If Y can not be cleaved and contains at least two points, then it is a chain of X 2 .

Proof. Let Y be a closed subset of X 2 that contains at least two points and can not be cleaved. The proof that Y is a chain has three steps.

Step 1:

Let Y ω k be any ring of X 2 , such that Y intersects both Y ω k and (Y ω k ) c . Then, Y ω k is contained in Y .
We prove this fact by contradiction. If Y ω k is not contained in Y , one can apply Lemma 2.3: the set A 1 = Y ∩ Y ω k may be contracted to a small euclidean ball contained in Int(T ω k ). Since (Y ) c is connected, one can then send this ball outside any bounded neighborhood U of T (0) , by an isotopy that fixes the closed set Y \ A 1 . The closed set A 2 = Y \ A 1 is not empty by assumption and it is possible to contract it in a small euclidean ball through an isotopy whose support is contained in U. This shows that Y can be cleaved. This is a contradiction and the fact is proved.

Step 2: There exists a word ω k such that Y is a union of some rings of the form Y ω k ,ℓ (with the same level k + 1). Let z 1 be any point in Y . By assumption, Y contains at least an other point z 2 . By construction of X 2 , there exists a ring Y ω k that contains z 1 and not z 2 . By the first step, Y ω k is included in Y . We hence proved that Y is a union of rings of X.

Let us note that, for any two rings A 1 and A 2 of X 2 , either they are disjoint or one is contained in the other. Consequently, Y is a disjoint union of rings that are maximal in Y for the inclusion.

Let A ⊂ Y be one of these maximal rings and let A ′ be the ring of X 2 which contains A and was built at the previous level. Hence, A ′ is not contained in Y . By the first step of this proof, one deduces that Y is contained in A ′ . Let us denote A ′ by T ω k . Repeating this argument with any maximal subring of Y , one deduces that Y is a union of rings of the form Y ω k ,i . This proves the second step.

Step 3: The set Y is a chain. From the second step, Y decomposes as an union of chains made of rings of the form Y ω k ,i . If Y is not a chain, it decomposes as chains C 1 , . . . , C r , that are not linked together. Hence, one can contract C 1 and C 2 ∪ • • • ∪ C r in two disjoint euclidean balls. This is a contradiction since Y can not be not cleaved. Hence, Y is a single chain.

End of the proof of Proposition 2.1. Let Y be a closed subset of X 2 that satisfies the properties of Proposition 2.1.

By the second property, Y decomposes as a union A 1 ∪ A 2 ∪ A 3 . By Lemma 2.4, all the sets A i are chains of X 2 . By the same lemma, the unions A i ∪ A j are chains as well, hence A 1 , A 2 and A 3 are unions of rings of the form Y ω k ,i contained in a same ring Y ω k . But since they are disjoint the only way that A 1 ∪ A 2 , A 2 ∪ A 3 and A 3 ∪ A 1 are all chains of X 2 is that their union is the ring Y ω k . This shows that Y is a ring and concludes the proof.

Rigidity of X 2 under homeomorphisms

We now end the proof of Theorem 1.3.

Let f be any homeomorphism from

R 3 onto itself such that f (X 2 ) ⊂ X 2 . From the Proposition 2.1, the image of any ring Y ω k by f is a ring Y τ j . Lemma 2.5. Let Y ω k be a ring of X 2 and Y τ j its image by f . The images of the subrings Y ω k ,i of level k + 1 of Y ω k are subrings Y τ j ,ℓ of level j + 1 of Y τ j .
Proof. Let us assume that it is not the case. There would exists a ring Y of level j + 1 in Y τ j which contains strictly the image of some subring Y ω k ,i of Y ω k . In other words, the preimage f -1 (Y ) is a proper subset of Y ω k but is larger than the subrings of level k + 1 of Y ω k . However, by Proposition 2.1, the set f -1 (Y ) is a ring of X 2 . This is a contradiction. By this lemma, the rings Y ω k and Y τ j decompose in the same number of subrings of next level so that n ω k = n τ j . By the construction, the lengths of the closed chains of tori in X 2 are all different. This implies that Y ω k = Y τ j . Since any point z ∈ X 2 is the intersection of some family of rings Y ω k (with k ր ∞), it is fixed by the map f . Hence, the restriction of f to X 2 is the identity. This concludes the proof of Theorem 1.3.

Proof of Proposition 1.4

Let X be a Cantor set on the plane (the proof for the Cantor set on the line is similar but simpler and will be left for the reader as an exercise). Let Σ = ∞ n=0 {0, 1} n be the space of all finite zero-one sequences.

Preliminary constructions

We start by the construction of the covering of X with a family of topological closed balls B ω n with ω n ∈ Σ, satisfying the following properties:

i) ∀ i∈{0,1} B ω n ,i ⊂ Int B ω n , ii) ∀ ω n ∈Σ ∀ i∈{0,1} B ω n ,0 ∩ B ω n ,1 = ∅, iii) ∀ n≥0 X ⊂ ω n ∈{0,1} n Int B ω n , iv) ∀ ω n ∈Σ B ω n ∩ X = ∅, v) lim n→∞ sup ω n ∈{0,1} n diam B ω n = 0.
We can easily do this construction for some special Cantor sets (like the usual middle one third Cantor set on a line, contained in the plane). As any two embeddings of Cantor sets in the plane are equivalent with respect to plane homeomorphisms (see [M], chap. 13, theorem 6, p. 93), we get the construction for X. As the distance from X to the boundary of any B ω n and the distance between boundaries of B ω n and B ω n i is non-zero, we may freely assume that the boundary of all balls

B ω n is C 1 . Let B (n) = ω n ∈{0,1} n B ω n
. By i), they form a decreasing sequence, by iii), iv) and v), we have X = n B (n) .

In the course of the proof, we will construct inductively another family of topological closed balls (even euclidean balls, but it is not really necessary) C ω n and define the contractive map F in such a way that F ((B ω n ) c ) = (C ω n ) c . The sets C ω n will satisfy the properties i), ii), iv) and v) above.

We define the sets C (n) similarly to B (n) : it is a decreasing family as well but its limit is only some Cantor subset of X.

Let us note that the construction of the family B ω n is precisely the point where the proof fails in dimension higher than two. In fact for the Cantor sets in R 3 for which such a construction is still possible, the rest of the proof follows with minor modifications. These Cantor sets can be easily seen as images of a Cantor set on the line R 1 ⊂ R 3 under homeomorphisms of R 3 .

A geometric lemma

We will use the following lemma: Lemma 2.6. Let A, A 0 , A 1 be closed topological balls with C 1 boundary, such that A 0 and A 1 are disjoint and contained in Int(A). Let r be a positive real number. Let z, z 0 , z 1 be three points such that z 0 , z 1 ∈ U = B(z, r).

Then, there exists a constant k(A, A 0 , A 1 , U, z 0 , z 1 ) with the following property: for any Lipschitz homeomorphism g : ∂A → ∂B(z, r) with Lipschitz constant L and for all sufficiently small r 0 , r 1 one can find a map h

: A \ Int(A 0 ∪ A 1 ) → B(z, r) \ Int(B(z 0 , r 0 ) ∪ B(z 1 , r 1 )) such that: a) h |∂A = g, b) h is a homeomorphism, h(∂A i ) = ∂B(z i , r i ),
c) h is Lipschitz and its Lipschitz constant is not greater than L ′ = L • k(A, A 0 , A 1 , U, z 0 , z 1 ), d) the Lipschitz constant of h |∂A i is not greater than

L ′′ = L • r i • k(A, A 0 , A 1 , U, z 0 , z 1 ).
Note that the constant k(A, A 0 , A 1 , U, z 0 , z 1 ) will not change if we rescale the triple {U, z 0 , z 1 } by a similitude or exchange z 0 with z 1 .

Proof. As all the possible triples (A, A 0 , A 1 ) are bi-Lipschitz equivalent, it is enough to prove the lemma for A = B(z, r), A 0 = B(z 0 , ρ), A 1 = B(z 1 , ρ) where ρ is the greatest such number that B(z 0 , 4ρ) and B(z 1 , 4ρ) are disjoint and contained in B(z, r). Similarly, we may assume that r = 1, z = (0, 0) (the rescaling changes L, L ′ and L ′′ by the same multiplicative constant) and that g is orientation preserving, ie. g(cos φ, sin φ) = (cos g 0 (φ), sin g 0 (φ))

for some orientation-preserving homeomorphism g 0 of S 1 . Note that now L ≥ 1. We assume r 0 , r 1 < ρ.

We define h as follows:

• on the annulus 1 ≥ |x| ≥ 1r by the formula h(a cos φ, a sin φ) = a cos 1a r φ + a + r -1 r g 0 (φ) , a sin 1a r φ + a + r -1 r g 0 (φ) .

• on the set B((0, 0), 1r) \ (B(z 0 , 2r) ∪ B(z 1 , 2r)) as the identity,

• on B(z i , 2r) by the formula h(z i + a • (cos φ, sin φ)) = z i + 2rr i r a + 2r i -2r • (cos φ, sin φ).

This map satisfies a) and b), its Lipschitz constant is not greater than L in the first part, 1 in the second one and 2 in the third one (this gives c)) and its Lipschitz constant on ∂A i is not greater than r i /r (we get d)).

Supported by Foundation for Polish Science and by KBN Grant 2 PO3A 034 25.

End of the constructions

We are now prepared to construct inductively the sets C ω n and the contractive homeomorphism F on each set Int(C (n) ) \ Int(C (n+1) ): we consider the Cantor set X, two points x and y belonging to X and the covering {B ω n } constructed in the claim above. Without weakening the assumptions, we may choose zeros and ones in our symbolic notation in such a way that {x} = B 0 n .

We begin with the definitions of C ∅ = C (0) and F on R 2 \ Int(C ∅ ). As the Cantor set is perfect, y is not an isolated point, hence we have a sequence (y ℓ ) in X converging to y. For any ℓ, the triple {B(y, 2|y ℓ -y|), y, y ℓ } is geometrically identical up to a similitude. We take any Lipschitz homeomorphism f 0 from B c ∅ onto B c (y, 2|y 0 -y|). We can construct a family of Lipschitz homeomorphisms

, where S ℓ is the similitude moving {B(y, 2|y 0 -y|), y, y 0 } onto {B(y, 2|y ℓ -y|), y, y ℓ }. The Lipschitz constants of those homeomorphisms are decreasing to zero together with |y ℓ -y|, hence for some ℓ they will be smaller than

Now we explain how to choose C 0 and C 1 and extend F on Int(C ∅ ) \ Int(C 0 ∪C 1 ). At this step, F is a contraction from the complement of Int(B ∅ ) onto the complement of Int(C ∅ ) and (by Lemma 2.6) it can be extended to an homeomorphism from the complement of Int(B 0 ∪ B 1 ) onto the complement of Int(B(y, r 0 ) ∪ B(y ℓ , r 1 )) for some small r 0 and r 1 . By Lemma 2.6 c) and by our choice of the Lipschitz constant of F in restriction to the boundary ∂B ∅ , the map F is a contraction.

Moreover by Lemma 2.6 d), for sufficiently small r 0 and r 1 the Lipschitz constant on ∂B 0 and ∂B 1 is arbitrarily small. We may thus choose a point from X very close to y and repeat the procedure, obtaining a contraction on Int(B 0 ) \ Int(B 00 ∪ B 01 ). Similarly choosing a point sufficiently close to y n lets us extend our function on Int(B 1 ) \ Int(B 10 ∪ B 11 ) and so on.

In this inductive procedure we will get a contracting map from the complement of X onto the complement of some Cantor set Y on the plane (the closure of the set of all the points we have chosen on all the stages of the construction). As the points chosen were always belonging to X, the set Y is included in X. We then extend the function F on X by continuity and we are done.


Sever Giovanni Manfredi
email: giovanni.manfredi@ipcms.u-strasbg.fr

Adrian Hirstoaga
email: hirstoaga@math.unistra.fr

Stéphane Devaux
email: stephane.devaux@ccfe.ac.uk

S Hirstoaga

Vlasov modelling of parallel transport in a tokamak scrape-off layer

come   L'archive ouverte pluridisciplinaire

Introduction

One of the main challenges for future tokamak operation, such as ITER, is constituted by the large heat load on the divertor plates. The divertor surfaces are constantly bombarded with high-energy neutral and charged particles and may thus see their lifetime considerably reduced [START_REF] Janeschitz | [END_REF]2]. The erosion also releases high-Z impurities, which migrate towards the bulk plasma and, due to radiation, deteriorate its confinement [START_REF] Stangeby | The plasma boundary of magnetic fusion devices[END_REF]. In order to keep within reasonable limits the erosion of different wall materials, it is important to estimate the plasma characteristics in the scrape-off layer (SOL), i.e. the region outside the last closed magnetic surface (separatrix). Several numerical and analytical studies have addressed the problem of the transition between a hot plasma and a material surface, both for unmagnetized and magnetized plasmas [START_REF] Chodura | [END_REF]5,6,7].

Most such studies focus on the static plasma-wall transition. However, transient events -such as edge-localized modes (ELMs) -routinely occur in the tokamak edge during H-mode confinement [8]. ELMs are plasma relaxations, presumably of MHD origin, which cause a sudden drop in density and temperature of the pedestal plasma, leading to a significant loss of the stored plasma energy (5 -20%) [9,[START_REF] Pitts | JET EFDA contributors[END_REF]. Once the ELM-driven plasma pulse has crossed the magnetic separatrix, it travels mainly parallel to the magnetic field lines and ends up hitting the divertor plate, after a delay typically of the order of a few hundred microseconds. Such violent events can pose a serious threat to the long-time resistance of the divertor materials.

The SOL is a complex region where a multitude of physical and chemical processes take place, occurring on rather disparate space and time scale. A realistic description of the SOL thus requires sophisticated modelling that takes into account, among others, plasma turbulence, collisions and ionisation, impurity transport, and plasma-surface interactions. The corresponding simulations [11,12] are computationally very costly, so that a thorough analysis of the various physical regimes is a demanding task.

On the other hand, a few semi-analytical models have been developed in order to capture the essential features of parallel transport in the SOL. In particular, an elegant and exactly-solvable model can be obtained by neglecting the Coulomb force in the parallel kinetic equation [13]. Together with an initial Maxwellian distribution, this surprisingly simple model reproduces with good accuracy some of the main features of an ELM signal, most notably its rapid rise (∼ 200µs) followed by a much slower decay (up to 3ms). The main drawback of this model is that it fails to respect quasineutrality, which is a direct consequence of the Coulomb attraction between the ions and the electrons. (Some form of 'weak' quasi-neutrality can be imposed a posteriori to partially correct this shortcoming [13]).

In the present paper, we adopt an intermediate stance: the fully nonlinear parallel dynamics is solved by means of a self-consistent kinetic model (Vlasov-Poisson equations), but we neglect collisions and effects due to plasma-surface interactions, such as secondary electron emission and recycling. The Vlasov equation is integrated numerically using an Eulerian method (meshing of the entire phase space), which guarantees a high resolution in velocity space, even in regions of low plasma density such as the sheaths [14]. In addition, thanks to a recently-developed 'asymptotic-preserving' numerical scheme [15,16], it is possible to lift numerical constraints on the time step and grid spacing, which are no longer limited by, respectively, the electron plasma period and Debye length.

The resulting numerical code is computationally manageable on a desktop computer (2-3 hours of CPU for a typical run, albeit with an effective parallel connection length smaller than its realistic value) and capable of reproducing with good accuracy the principal features of parallel plasma transport during ELMs. The remainder of this paper will be devoted to the description of the model and the analysis of the numerical results for some experimentally-motivated scenarios.

The Vlasov-Poisson approach

In our approach, we assume that the charged particles (or rather, their guiding centres) travel along the magnetic field lines, but not across them. Thus, we can adopt a onedimensional geometry along the parallel direction, here denoted x, with corresponding parallel velocity v x . In the perpendicular plane, the distribution function remains Maxwellian at all times, so that the distribution in the four-dimensional phase space (x, v) reads as:

F j (x, v, t) = f j (x, v x , t)M j (v ⊥ ), (1)
where

M j (v ⊥ ) = (1/2πv 2 T j ) exp(-v 2 ⊥ /2v 2 T j
), v T j = T j /m j is the thermal speed, and the label j = i, e stands for ions and electrons respectively. The temperature is supposed to be the same for both species, T i = T e = T ELM .

Under these assumptions, the ion and electron evolutions are described by the onedimensional Vlasov equations in the parallel phase space (x, v x ):

∂f j ∂t + v x ∂f j ∂x - q j m j ∂φ ∂x ∂f j ∂v x = S j (x, v x , t), (2)
where q j = ±e (we consider singly charged ions). The collisionless approximation is reasonable for the initial phase of high-energy transients such as ELMs, for which the thermal mean-free-path exceeds the parallel connection length L , defined as the typical distance between the outer midplane and the divertor outer target. In the JET tokamak, the parallel connection length is roughly L ∼ 30 m. The source terms S j describe the plasma pulse generated by the ELMs, and can be written as S j (x, v x , t) = f ELMj (x, v x )g(t), where

f ELMj (x, v x ) = n ELM exp - x 2 2σ 2 exp(-v 2 x /2v 2 T j ) √ 2πv T j , (3)
and g(t) models the pulse temporal profile. Thus, the source is a Maxwellian distribution, whose spatial profile is given by a Gaussian with width σ ∼ 0.1L and peak density n ELM . Finally, the electrostatic potential φ(x, t) obeys the one-dimensional Poisson equation

∂ 2 φ ∂x 2 = - e ǫ 0 (n i -n e ) . (4)
The above equations are solved on an interval x ∈ [-L , L ], where x = ±L represent the locations of the divertor plates. The plates are supposed to be perfectly absorbing surfaces (i.e. the incoming flux is zero) and are kept at constant electric potential, φ(±L ) = 0. Besides being relevant to the ELM dynamics in a tokamak, the above model describes the general scenario of an initially neutral plasma created in a confined region delimited by two grounded plates: the plasma first undergoes an expansion into vacuum and then reaches the plates. Here, we are mainly interested in the plasma features on the plates, particularly the particles and energy fluxes, defined respectively as (we omit the species index for clarity):

Γ(t) = v x F (±L , v, t)dv, (5)
Q(t) = 1 2 m(v 2 x + v 2 ⊥ )v x F (±L , v, t)dv = Q + Q ⊥ . (6)
It can be readily checked that Q ⊥ = Γ T ELM . We use typical parameters for Type-I ELMs [13]: plasma density n ELM = 5 × 10 19 m -3 , temperature T ELM = 1.5keV, and parallel connection length L = 30m. With these values (and considering hydrogen ions, m i = 1836m e ) the ion traversal time is τ i = L /v T i ∼ 80µs and the Debye length λ D ∼ 40µm. In the simulation results, the parallel energy fluxes will be expressed in terms of Q ref = n ELM v T i T ELM , which, for the above set of parameters, is

Q ref = 4.5GWm -2 .
The relevant physical regime will be determined by two dimensionless parameters, namely, the ion-to-electron mass ratio and the ratio of the Debye length to the parallel connection length, λ ≡ λ D /L ∼ 10 -6 . The latter is particularly important as it determines the scale length over which quasi-neutrality can be violated. This can be easily seen by rewriting Poisson's equation in dimensionless units, normalizing space to L , the densities to n ELM , and the electric potential to T ELM /e:

λ 2 ∂ 2 φ ∂x 2 = n e -n i . (7)
When λ → 0, the only way to satisfy Poisson's equation is for the plasma to be quasineutral (n e -n i → 0). The electric potential does not vanish, but takes exactly the value required to enforce the quasi-neutrality constraint. In that case, Poisson's equation cannot be used to determine the potential, as Eq. ( 7) becomes singular. When λ is finite but very small, as is the case for our problem, quasi-neutrality is satisfied over distances much larger than the Debye length, but breaks down at finer scales (a typical instance are the sheaths that form at the plasma-surface boundary, which have a thickness of a few Debye lengths). Thus, an increasingly fine spatial resolution is required as λ → 0. In addition, since the plasma frequency scales as λ -1 , the time step must also be reduced in the same proportions. Indeed, it has been proven that a numerical instability occurs if the adopted time step is larger than the inverse plasma frequency [15].

For the above reasons, it would be useful to replace the standard Poisson's equation (7) with an equivalent equation that can provide φ when λ → 0. This can actually be achieved rather easily [15,16] by using the mass and momentum conservation equations derived from the Vlasov equations ( 2). The resulting 'reformulated' Poisson's equation, written in appropriate dimensionless units, reads as:

λ 2 ∂ 2 ∂t 2 ∂ 2 φ ∂x 2 + ∂ ∂x (ε n i + n e ) ∂φ ∂x = - ∂ 2 ∂x 2 (R i -R e ) , (8)
where ε = m e /m i and R j = f j v 2 dv (see the Appendix A for details on its derivation). Equation ( 8) is completely equivalent to the standard Poisson equation (provided the latter and its time derivative are satisfied at t = 0), but does not become singular in the asymptotic regime λ → 0. For this reason, the above approach has been termed 'asymptotic-preserving' method.

Thanks to the asymptotic-preserving scheme, it is possible to employ a grid spacing that exceeds the Debye length and, more importantly, a time step larger than the inverse plasma frequency. In the forthcoming simulations, we have used a grid spacing ∆x = 2λ D and a time step up to ∆t = 4ω pe . Note that here λ D and ω pe refer to their initial values computed with the peak density n ELM and temperature T ELM . These values, in the normalisation of Eq. ( 7), are equal to λ and λ -1 , respectively. Some details on the numerical techniques used to solve the model equations are provided in Appendix B.

Free-streaming model

Fundamenski et al. [13] have developed a simple analytical model for parallel transport, which nevertheless captures the main features of an ELM signal, most notably its rapid rise and slower decay. The model completely neglects Coulomb interactions, so that the Vlasov equations are reduced to the following free-streaming equations (in the rest of this section we omit the species index for simplicity):

∂f ∂t + v x ∂f ∂x = S(x, v x , t). (9)
This equation possesses the exact solution

f (x, v x , t) = t 0 f ELM (x -v x t + v x t ′ , v)g(t ′ )dt ′ . ( 10
)
We further assume that the source is instantaneous in time, g(t) = δ(t), where δ is the Dirac delta function, and that it is localized in space,

n ELM (x) ≡ ∞ -∞ f ELM (x, v)dv = √ 2πσδ(x)
. With these hypotheses, it is possible to express analytically all relevant quantities at an arbitrary location x.

For the particles density, we obtain:

N (x, t) ≡ F (x, v, t)dv = σ v T n ELM t exp - x 2 2t 2 v 2 T , (11)
whereas the particles flux yields: Γ = (x/t)N . Further, we have for the kinetic energy:

E(x, t) ≡ 1 N F (x, v, t) m 2 v 2 dv = E + E ⊥ = m 2 x 2 t 2 + T ELM (12)
and the energy flux:

Q(x, t) = Q + Q ⊥ = m 2 x 2 t 2 Γ + ΓT ELM . ( 13
)
Using these expressions, we can define the parallel sheath transmission coefficient: γ = Q /(ΓE ), which is exactly equal to unity within this model.

The above formulae will be used as a benchmark for the numerical results presented in the next section.

Simulation results

We now present the numerical results obtained by solving the full Vlasov-Poisson equations, i.e. including the effect of the electric field. We shall use two types of temporal profiles for the sources in the Vlasov equations ( 2): an instantaneous source g(t) = δ(t), and a distributed-in-time source with

g(t) = Ct 2 exp - (t -t 0 ) 2 2σ 2 t , (14)
where the normalization constant C is chosen so that ∞ 0 g(t)dt = 1. Provided that σ t ≪ t 0 , the above expression peaks at t = (t 0 /2)(1 + 1 + 8σ 2 t /t 2 0 ) and has the advantage of vanishing at t = 0.

In all cases, except where otherwise stated, we use the following parameters: m i = 1836m e , σ i = 0.1L , and λ = 10 -3 . This value of λ is much larger than the realistic value for ELM plasmas. However, we have checked that the evolution of relevant quantities (such as the particles and energy fluxes) are virtually unchanged when smaller values of λ are taken (an exception is an early peak in the electron energy flux, which disappears in the asymptotic limit λ → 0; see later for further details). Thus, it appears that for λ = 10 -3 we have already reached the asymptotic regime and the results can be confidently extrapolated to the realistic value of λ.

Instantaneous source

First of all, we tested our numerical code by comparing the analytical results from the free-streaming model with the numerical solution of the Vlasov equation without electric field (9). As expected, the analytical and numerical fluxes (not shown here) were virtually identical.

Next, we turn to the full solution of the Vlasov-Poisson equations including the electric field. Figure 1 shows the evolution of the electron and ion particles and energy fluxes on the divertor plates (since this case is symmetric there is no difference between the inner and the outer plate), for λ = 10 -3 and λ = 2.5 × 10 -4 . The standard pattern (sharp rise, followed by a slower decay) is recovered in the simulations.

Very little difference in the fluxes is observed between the λ = 10 -3 and the λ = 2.5 × 10 -4 cases. Therefore, it seems that the asymptotic regime (λ → 0) has already been reached, so that the results for the fluxes should remain unchanged for a realistic value of λ, which is of the order 10 -6 .

The ion and electron particles fluxes are virtually identical, which is a consequence of quasi-neutrality. In contrast, the electron energy flux is appreciably smaller than that of the ions. This is in line with measurements of the heat load on a tokamak divertor plates, where the electron load is observed to be roughly one third of the ion load [17,[START_REF] Eich | JET EFDA contributors, and ASDEX Upgrade team[END_REF]. The free-streaming model (dotted line in Fig. 1) is fairly accurate for the ion fluxes, especially over long times.

The above dynamics occurs on the timescale of the ionic transit time τ i , as is expected from quasi-neutrality. However, an early burst of electrons is observed on the electronic transit time scale, as can be seen from Fig. 2. This burst is due to a bunch of electrons escaping the attraction of the ions and reaching the plates before quasineutrality can be established. Indeed, the height of this peak decreases with decreasing λ, i.e. when quasi-neutrality is stronger. From the free-streaming result [Eq. ( 13)], the electron energy flux should peak at t peak = L/2v T e ∼ 0.5τ e . The observed peak is at roughly one half of this value, signalling that the peak is caused by suprathermal electrons travelling with an average velocity equal to twice the thermal speed.

It would be interesting to deduce a scaling law for the deposited electron energy as a function of λ. This can be done by integrating the electron fluxes on Fig. 2 between t = 0 and, say, t = τ e . The electron kinetic energy deposited on the divertor plate (x = L ) is then defined as E elec = τ e 0 Q e (t)dt and is shown in Fig. 3 on a log-log plot. It is clear that E elec depends quadratically on λ. In the present case, we obtain a scaling law of the type: E elec /T ELM ≈ 15λ 2 , which allows us to extrapolate the result to realistic values of λ.

The mechanism responsible for the above effect is best understood by inspecting the parallel phase space (x, v x ), displayed in Fig. 4, where we show the contour plots of the electron distribution function at various times. The electron velocity distribution on the right boundary is shown on Fig. 5 at the same instants. Very early, a small fraction (about 0.2%) of energetic electrons is ejected from the bulk plasma and hits the divertor plate. This leaves behind a net positive charge that traps the remaining electrons in a potential well, preventing them to reach the divertor plates. The potential well expands slowly following the motion of the ions (on a timescale ∼ τ i ) and finally reaches the boundary, thus releasing the ions' and the (remaining) electrons' energy onto the divertor plate. On the timescale of the ion dynamics, the velocity distributions on the wall become approximately half-Maxwellian for the electrons and Maxwellian for the ions [19], as can be seen from Fig. 6.

The densities and temperatures of the two species of particles on the plate are depicted in Fig. 7. The temperatures are particularly interesting: whereas the electron temperature converges to a constant (about one thousandth of the initial value T ELM ), the ion temperature decreases with time as t -2 . This type of decrease is readily obtained from the free-streaming model and was anticipated in earlier studies [20]. (One should keep in mind that here we are talking about the parallel temperatures; the perpendicular temperatures remain constant for all times, since no coupling between the parallel and perpendicular motion is allowed in our model).

The parallel sheath transmission coefficients are shown in Fig. 8. Whereas γ i quickly reaches unity, as expected from the free-streaming model, γ e is twice as large. This is due to the fact that the electronic distribution on the divertor plate is a half-Maxwellian (see Fig. 6), so that the parallel kinetic energy E is twice as small compared to a full Maxwellian. Finally, we show the time evolution of some integrated quantities, particularly the kinetic and potential energies (Fig. 9). As long as the main plasma population does not reaches the plates, it expands freely into a vacuum. During this phase, lasting up to ∼ 0.2τ i , some kinetic energy is transferred from the electrons to the ions.

This transfer occurs because of the fast electrons that move rapidly forwards leaving the ions behind, as was shown in Figs. 45. This creates a charge imbalance -and therefore an electric field -at the plasma-vacuum front, which is revealed by the peak of the potential energy in Fig. 9. The electric field slows down the electrons and accelerates the ions, thus effectively transferring some kinetic energy from the former to the latter. Eventually, the ions catch up with the electrons, the plasma becomes locally neutral, and the potential energy goes to zero. Over longer times (t > 0.2τ i ), the charged particles start being lost on the divertor plates and consequently both kinetic energies decrease.

Time-distributed source

In actual ELM events, the plasma pulse is not instantaneous, but has a finite duration, which is of the order of a few ion transit times. In order to simulate this type of scenario, we use a pulse temporal profile such as that of Eq. ( 14), with σ t = 0.7τ i and t 0 = 1.4τ i [the plot of g(t), given by the dotted line on Fig. 10, shows that the pulse peaks at about 1.9τ i ].

The evolution of the particles and energy fluxes is shown in Fig. 10. The fluxes reach a maximum at roughly one ion transit time after the source peak. The ion energy flux is about 30% larger than that of the electrons, whereas it was almost four times larger in the case of an instantaneous source (Fig. 1). It is useful to compare these results to a free-streaming case, for which the electric field was neglected (Fig. 11). In the free-streaming case, the electron fluxes peak earlier, basically in phase with the source (at closer inspection, the electron fluxes peak around 2-3 electron transit times after the peak of the source). Also, a larger fraction of the energy flux is carried by the electrons, at the expense of the ions' energy flux. Clearly, the electric field has the effect of transferring some energy from the electrons to the ions.

By inspecting the electron energy flux for short times (Fig. 12), we also observe a small peak occurring on the electronic timescale, as we had found for the case of an instantaneous source. Now the peak appears later, since the plasma needs to build up progressively, but still occurs much earlier than the peak of the ELM pulse (see inset of Fig. 12).

The velocity profiles on the divertor plate (Fig. 13) are similar to those obtained in the case of an instantaneous source, i.e. half-Maxwellian for the electrons and roughly Maxwellian for the ions.

Finally, the parallel sheath transmission coefficients are shown in Fig. 14, revealing some interesting behaviour. The ion coefficient remains close to unity for the entire duration of the run, signalling that the Vlasov-Poisson result stays close to the free-streaming solution. The electronic coefficient displays a long transient, which corresponds to the time-window when the source is active, then quickly relaxes to its asymptotic value (γ e = 2) after the source is extinguished.

Conclusion

Parallel transport plays a major role in the dynamics of ELM events in a tokamak. Even relatively simple models that ignore Coulomb interactions [13] manage to reproduce with fair accuracy the main features of an ELM event, particularly the energy load on the divertor plates.

Stimulated by this observation, we constructed a minimal kinetic model of parallel transport during ELMs, based on the Vlasov-Poisson equations in the direction parallel to the magnetic field (in the perpendicular plane, the plasma is supposed to be Maxwellian). The results of numerical simulations performed with this model were compared systematically to the free-streaming solutions [13].

Several interesting conclusions can be drawn from our study:

(i) Even for moderate connection lengths (we used L = 1000λ D throughout this work) the particles and energy fluxes are correctly described, i.e. by decreasing the ratio λ = λ D /L the fluxes are virtually unchanged;

(ii) Whereas the ion fluxes are relatively well reproduced by the free-streaming model, the electron fluxes are not. This is because the free-streaming model does not properly include quasi-neutrality;

(iii) On the electron transit time scale, an early burst can be observed on the divertor plates, which corresponds to suprathermal electrons escaping the ions' attraction due to their large kinetic energy. The remaining electrons are trapped in the potential well created by the ions (which now outnumber the electrons) and evolve in unison with them. On the ion timescale, the plasma is everywhere neutral, except for a positive net charge in front of the divertor plates (sheaths);

(iv) The long-time velocity distributions on the plates are Maxwellian for the ions and half-Maxwellian for the electrons;

(v) The parallel sheath transmission coefficients tend rapidly to 1 for the ions and 2 for the electrons;

(vi) Most of the above results hold equally for an ELM pulse with a realistic duration t ELM ∼ τ i .

The present model could be further improved in several ways. First, a pre-ELM environment (low temperatures and densities) should be computed self-consistently and used as a starting background prior to the ELM pulse. Second, even though collisions are not dominant due to the high temperature of the ELM pulse, they may play a role over long times, so that a suitable collision term should be added to the Vlasov equations.

Finally, the perpendicular dynamics should be treated with more sophistication, by taking into account that the magnetic field decreases with increasing distance from the tokamak axis. This introduces a 'mirror' force that acts differently on ions with different magnetic moments (this effect should be negligible for the electrons) and couples the parallel and perpendicular dynamics. A sufficiently accurate description could be achieved by assuming that the magnetic moment is an exact invariant and then using an approach such as that of Ref. [21]. The above issues are currently under investigation.

As mentioned in Section 2, in order to compute the potential φ leading to a quasineutral regime (the limit λ → 0) we need a reformulation of the Poisson's equation that does not become singular in this limit. Such an equation was already derived in the case of the Euler-Poisson [15] and the Vlasov-Poisson systems [22,16]. For the sake of completeness, we show how the reformulated Vlasov-Poisson model is obtained in our case.

We start with the following scaled Vlasov-Poisson equations

∂f i ∂t + v x ∂f i ∂x -ε ∂φ ∂x ∂f i ∂v x = S i (x, v x , t) (A.
= x L v x = v x v Te t = v Te L t n j = n j n ELM φ = φ T e /e (A .4)
(the overcarets, which will be omitted in the sequel, stand for the dimensionless variables). By taking the first two moments of the difference between Eq. (A.1) and Eq. (A.2), we obtain

∂ t (n i -n e ) + ∂ x (J i -J e ) = 0 (A.5) ∂ t (J i -J e ) + ∂ x (R i -R e ) + (ε n i + n e ) ∂ x φ = 0, (A.6)
where J j = f j (x, v x , t) v x dv x and R j = f j (x, v x , t) v 2 x dv x . Now we subtract the spatial derivative of Eq. (A.6) from the time derivative of Eq. (A.5):

∂ 2 ∂t 2 (n i -n e ) - ∂ ∂x (ε n i + n e ) ∂φ ∂x = ∂ 2 ∂x 2 (R i -R e ) . (A.7)
Finally, using Poisson's equation we express ∂ 2 tt (n i -n e ) in terms of φ and thus recover the reformulated Poisson equation (8).

It can further be shown [16,22] that the reformulated Poisson equation (8) (coupled to the Vlasov equation) is rigorously equivalent to the standard Poisson equation (7) provided that φ satisfy the following conditions at the initial time: (λ 2 ∂ 2 xx φ+n i -n e ) t=0 = 0 and ∂ ∂t (λ 2 ∂ 2 xx φ + n i -n e ) t=0 = 0.

Appendix B. Numerical methods

The Vlasov equation ( 2) is solved using an Eulerian method based on a uniform meshing of the parallel phase space (x, v x ). For the time-stepping, a second-order splitting scheme is used [14], which solves alternatively the advection in real space, ∂ t f + v x ∂ x f = 0, and the advection in velocity space, ∂ t f +(qE/m)∂ vx f = 0. Each advection step is performed using a third-order positive flux-conservative (PFC) method, with a slope corrector that prevents the distribution function from becoming negative [23,24]. A typical simulation requires 1000 points both in real space and in velocity space.

As to the reformulated Poisson's equation, it was shown that an explicit timestepping scheme is only stable if ∆t ≤ ω -1 pe ∼ λ [15]. In order to lift this numerical constraint, the reformulated Poisson equation must be solved through an implicit scheme. More precisely, we first discretize Eq. ( 8) as follows: ) .

- ∂ ∂x λ 2 ∂ x φ m+1 -2∂ x φ m + ∂ x φ m-1 ∆t 2 + (εn m i + n m e ) ∂ x φ m+1 = ∂ 2 ∂x 2 (R m i -R m e ).
Next, we make use of a discrete form of Eq. which allows us to compute φ m+1 from known quantities at time t m [16]. A finite difference method is then used to approximate the spatial derivatives in the equation.

Although the above time-stepping strategy is only first order in time, the method could be easily extended to second order accuracy. Some work along this direction was presented in Ref. [15].

Finally, we remark that: (i) when ∆t → 0 for fixed λ, Eq. (B.1) reduces to a discrete form of the standard Poisson's equation; (ii) when λ → 0 for fixed ∆t, the scheme leads to a consistent (non-singular) discretization of the λ → 0 limit of Eq.

Figure 1 .

 1 Figure 1. Particles (left panel) and energy (right panel) fluxes as a function of time normalized to the ion transit time. Solid lines refer to the ions, dashed lines to the electrons. The dotted lines represent the ion fluxes obtained from the free-streaming model. Top panels: λ = 10 -3 ; bottom panels: λ = 2.5 × 10 -4 .

Figure 2 .

 2 Figure 2. Electron energy flux on the divertor plate for short times, for three values of the normalized Debye length: λ = 10 -3 (solid line), λ = 5 × 10 -4 (dashed line), and λ = 2.5 × 10 -4 (dotted line). Note that here time is normalized to the electron transit time τ e = (m e /m i )τ i .

Figure 3 .

 3 Figure 3. Diamonds: Electron kinetic energy deposited on the right plate (x = L ) as a function of λ. The straight solid ine represents the curve E elec /T ELM = 15λ 2 .

Figure 4 .

 4 Figure 4. Contour plots of the electron distribution function in the parallel phase space at different times, measured in units of the electron transit time τ e = (m e /m i )τ i . Only contour levels for which f e ≤ 10 -4 are shown.

Figure 5 .

 5 Figure 5. Electron parallel velocity distribution at x = L at different times, measured in units of the electron transit time τ e = (m e /m i )τ i .

Figure 6 .

 6 Figure 6. Electron (left panel) and ion (right panel) velocity distribution at x = L at the end of the run (t = 3.73τ i ).

Figure 7 .

 7 Figure 7. Time evolution of densities (left panel) and parallel temperatures (right panel) at x = L . The temperatures are normalized to their initial value T ELM . Solid line: ions; dashed line: electrons. The dotted line in the left panel represents the theoretical prediction of Eq. (11) for the ions. The inset in the right panel shows the ion temperature multiplied by t 2 .

Figure 8 .

 8 Figure 8. Parallel sheath transmission coefficients for the electrons (diamonds) and the ions (stars).

Figure 9 .

 9 Figure 9. Time evolution of the total ion kinetic energy (dotted line), electron kinetic energy (dashed line), potential energy (×1000) (thin solid line), and total energy (thick solid line). All quantities are normalized to the initial total energy.

Figure 10 .

 10 Figure 10. Particles (left panel) and energy (right panel) fluxes as a function of time, for a case with time-distributed source. Solid lines refer to the ions, dashed lines to the electrons. The dotted lines represent the source temporal profile g(t), in arbitrary units.

Figure 11 .

 11 Figure 11. Same as Fig.10for free-streaming particles (no electric field).

Figure 12 .

 12 Figure 12. Electron energy flux on the divertor plate for short times (in units of the electron transit time). The inset shows the temporal profile of the source g(t).

Figure 13 .

 13 Figure 13. Electron (left panel) and ion (right panel) velocity distribution at x = L at the end of the run (t = 9.34τ i ), for a case with time-distributed source.

Figure 14 .

 14 Figure 14. Parallel sheath transmission coefficients for the ions (solid line) and the electrons (dashed line), for a case with time-distributed source. The dotted lines represent the source temporal profile g(t), in arbitrary units.

  Then, using Poisson's equation, we replace ∂ 2 xx φ m and ∂ 2 xx φ m-1 by density terms at the corresponding time steps:

  (A.5) to replace (n m i -n m e ) -(n m-1 i -n m-1 e ) with -∆t (∂ x J m i -∂ x J m e ).We thus obtain the discrete (in the time variable) version of the reformulated Poisson's equation-∂ ∂x λ 2 + ∆t 2 (εn m i + n m e ) ∂ x φ m+1 = ∆t 2 ∂ 2 ∂x 2 (R m i -R m e ) + (n m i -n m e ) -∆t (∂ x J m i -∂ x J m e ), (B.1)

Acknowledgments

We thank W. Fundamenski for several useful discussions.

Appendix A. Reformulated Poisson's equation