Samy Blusseau

Santiago Velasco-Forero

Jesús Angulo

Isabelle Bloch

Morphological adjunctions represented by matrices in max-plus algebra for signal and image processing

Keywords: Morphological operators, Max-plus algebra, Graph theory

In discrete signal and image processing, many dilations and erosions can be written as the max-plus and min-plus product of a matrix on a vector. Previous studies considered operators on symmetrical, unbounded complete lattices, such as Cartesian powers of the completed real line. This paper focuses on adjunctions on closed hypercubes, which are the complete lattices used in practice to represent digital signals and images. We show that this constrains the representing matrices to be doubly-0-astic and we characterise the adjunctions that can be represented by them. A graph interpretation of the defined operators naturally arises from the adjacency relationship encoded by the matrices, as well as a max-plus spectral interpretation.

Introduction

Like linear filters can be represented by matrices in discrete image and signal processing, many morphological dilations and erosions can be seen as applying a matrix product to a vector, but in the minimax algebra. This is in particular the case for those defined with structuring functions, either flat or not, local or non-local [START_REF] Salembier | Study on nonlocal morphological operators[END_REF][START_REF] Velasco-Forero | On nonlocal mathematical morphology[END_REF], translation invariant or spatially variant [START_REF] Bouaynaya | Theoretical foundations of spatially-variant mathematical morphology Part I: Binary images[END_REF][START_REF] Debayle | General adaptive neighborhood image processing: Part I: Introduction and theoretical aspects[END_REF][START_REF] Lerallut | Image filtering using morphological amoebas[END_REF][START_REF] Verdú-Monedero | Anisotropic morphological filters with spatially-variant structuring elements based on image-dependent gradient fields[END_REF]. They are commonly known to be the vertical-shift-invariant dilations and erosions [START_REF] Heijmans | The algebraic basis of mathematical morphology I. dilations and erosions[END_REF][START_REF] Maragos | Chapter Two -Representations for morphological image operators and analogies with linear operators[END_REF]. While the matrix point of view is not the most appropriate for the implementation of these operators, especially translation-invariant ones, it is a valuable insight for their theoretical understanding. In particular, it can help predict and control complex behaviours such as those of iterated operators based on adjunctions with non-flat, spatially variant and input-adapted structuring functions [START_REF] Blusseau | Tropical and morphological operators for signal processing on graphs[END_REF][START_REF] Blusseau | Adaptive Anisotropic Morphological Filtering Based on Co-Circularity of Local Orientations[END_REF]. Indeed, it is a flexible and general framework which embraces a very broad part of morphological literature, and it is supported by the rich theory of Minimax algebra [START_REF] Akian | Max-plus algebra[END_REF][START_REF] Cuninghame-Green | Minimax algebra[END_REF].

In the abundant literature on spatially-variant morphological image processing, only a few approaches explicitly used the matrix formulation [START_REF] Blusseau | Tropical and morphological operators for signal processing on graphs[END_REF][START_REF] Blusseau | Adaptive Anisotropic Morphological Filtering Based on Co-Circularity of Local Orientations[END_REF][START_REF] Salembier | Study on nonlocal morphological operators[END_REF][START_REF] Velasco-Forero | On nonlocal mathematical morphology[END_REF], whereas most contributions were limited to flat structuring elements and focused on the local effects of the adaptive strategy. On the theoretical side, the representation of morphological adjunctions by matrices was studied in a setting that does not directly apply to digital signal and image processing, as the co-domain is usually an unbounded lattice, stable under any vertical translation [START_REF] Maragos | Chapter Two -Representations for morphological image operators and analogies with linear operators[END_REF]. Although a method was proposed to convert these adjunctions to new ones on bounded lattices [START_REF] Ronse | Why mathematical morphology needs complete lattices[END_REF], it is not practical and does not allow for the interpretations that are exposed here.

In the present paper, we focus on complete lattices of the type [a, b] n , where a and b represent the minimal and maximal possible signal values (typically a = 0 and b = 255 for 8-bits images), and n is an integer representing the size of the signal (typically, the number of pixels of an image, reshaped as a column vector). This is a theoretical contribution that can be viewed as a companion paper to previous studies where this framework has been successfully applied to adaptive anisotropic filtering [START_REF] Blusseau | Tropical and morphological operators for signal processing on graphs[END_REF][START_REF] Blusseau | Adaptive Anisotropic Morphological Filtering Based on Co-Circularity of Local Orientations[END_REF] 3 . In Section 2 we introduce the matrix-based morphological setting and prove simple but fundamental results: in particular, we characterise the adjunctions that can be represented by matrices and show that these matrices need to be doubly-0-astic. By viewing matrices as encoding adjacency, we provide in Section 3 a graph interpretation of iterated operators and their associated granulometries. In Section 4 we draw a link between these operators and some results on the spectrum of matrices in the max-plus algebra, before concluding in Section 5.

Matrix-based morphological adjunctions

Notations

In this paper matrices will be denoted by capital letters, such as W , and their i-th row and j-th column coefficients by corresponding indexed lowercase letters w ij . Similarly, vectors are written as boldface lowercase letters, such as x, and their i-th component as x i . Let 0 ≤ a < b ∈ R + be two non-negative real numbers, n ∈ N * a positive integer. The set {1, . . . , n} will be denoted by 1, n . Let L = ([a, b] n , ≤) be the complete lattice equipped with the usual product partial ordering (Pareto ordering): x ≤ y ⇐⇒ x i ≤ y i , ∀i ∈ 1, n . The supremum and infimum on L are induced by the Pareto ordering: for a family (x (k) ) k∈K of L, k∈K x (k) is the vector y defined by

y i = k∈K x (k)
i , where K is any index set. Therefore a = (a, . . . , a) T and b = (b, . . . , b) T are respectively the smallest and largest elements in L. For x ∈ L, we note x c = bx + a, and for any i ∈ {1, . . . , n}, e (i) is the "impulse" vector in L such that e

(i) i = b and e (i) j = a for j = i.
We note R max = R ∪ {-∞}, R min = R ∪ {+∞} and M n the set of n × n square matrices with coefficients in R max . Like (R max , ∨, +), (M n , ∨, ⊗) is an idempotent semiring, with the addition ∨ and product ⊗ defined as follows. For A, B ∈ M n , A ⊗ B and A ∨ B are the n × n matrices defined respectively by

(A ⊗ B) ij = n k=1 a ik + b kj and (A ∨ B) ij = a ij ∨ b ij = max(a ij , b ij ), for 1 ≤ i, j ≤ n. Similarly, for x ∈ R n max , A ⊗ x is the vector such that (A ⊗ x) i = n j=1 a ij + x j .
Note that ∨ and ⊗ are associative and ⊗ is distributive over ∨. Finally, the product of a scalar λ ∈ R max by a vector x ∈ R n max is λ ⊗ x = λ + x, the vector in R max such that (λ ⊗ x) i = λ + x i . In [START_REF] Cuninghame-Green | Minimax algebra[END_REF] and [START_REF] Velasco-Forero | On nonlocal mathematical morphology[END_REF], special subsets of M n are introduced, that we will show to be essential to represent morphological adjunctions on L.

Definition 1 (0-asticity [START_REF] Cuninghame-Green | Minimax algebra[END_REF]) A matrix W ∈ M n is said row-0-astic if for any 1 ≤ i ≤ n, n j=1 w ij = 0. Similarly, it is said column-0-astic if the supremum of each column is 0, and doubly-0-astic if the matrix is both row-0astic and column-0-astic. Finally, W is simply said 0-astic if 1≤i,j≤n w ij = 0.

A special kind of doubly-0-astic matrices are those with zeros on the diagonal and non-positive coefficients elsewhere.

Definition 2 (CMW matrices [START_REF] Velasco-Forero | On nonlocal mathematical morphology[END_REF]) A matrix W ∈ M n is called a Conservative Morphological Weights (CMW) matrix if ∀i, j ∈ 1, n , w ij ≤ 0 and w ii = 0.

We now introduce the morphological framework on L, based on the max-plus algebra product between matrices and vectors.

Dilations

For W ∈ M n , we consider the function

δ W from L to R n max such that ∀x ∈ L, δ W (x) = W ⊗ x =   1≤j≤n {w ij + x j }   1≤i≤n . (1)
In the processing of digital data such as images we usually want the input to be comparable with the output. Hence, we will constrain W such that δ W (L) ⊆ L. This has the following consequences:

δ W (b) ≤ b ⇒ ∀i ∈ 1, n , b + ( n j=1 w ij ) ≤ b ⇒ ∀i ∈ 1, n , n j=1 w ij ≤ 0 since b > -∞. Similarly, δ W (a) ≥ a ⇒ ∀i ∈ 1, n , n j=1 w ij ≥ 0.
Hence a necessary condition to have δ W (L) ⊆ L is that W be row-0-astic (Def. 1). Conversely, the row-0-asticity for W implies that δ W (a) = a and δ W (b) = b, and therefore that δ W (L) ⊆ L by increasingness of δ W . This leads to the following result.

Proposition 1 Let W ∈ M n and δ W be the function defined by [START_REF] Akian | Max-plus algebra[END_REF]. Then δ W is a dilation mapping L to L if and only if W is row-0-astic.

Proof. If δ W is a dilation mapping L to L, then δ W (L) ⊆ L which, as we showed, implies that W is row-0-astic. Conversely, we saw that a row-0-astic W implies δ W (L) ⊆ L. Therefore, we only have to verify that δ W is a dilation, or equivalently that it commutes with the supremum. This is straightforward from the definition of W ⊗ x. ⊓ ⊔

Erosions and adjunctions

Now we suppose that W ∈ M n is row-0-astic, hence δ W is a dilation from L to L, and we are interested in its adjoint erosion α W defined for any y ∈ L by α W (y) = E y where E y = {x ∈ L, δ W (x) ≤ y}. Let us denote by ε W the function from L to R n min such that for any y ∈ L

ε W (y) = δ W T (y c ) c = W T ⊗ y c c =   1≤j≤n {y j -w ji }   1≤i≤n . (2)
Then we can check that ∀y ∈ L, α

W (y) = ε W (y) ∧ b. Indeed, from (1) we see that for any x, y ∈ L, δ W (x) ≤ y ⇐⇒ x ≤ ε W (y). Therefore, since δ W (a) = a ≤ y we get a ≤ ε W (y), which implies ε W (y) ∧ b ∈ L; furthermore, ε W (y) ∧ b ≤ ε W (y) so ε W (y) ∧ b ∈ E y ; finally, as both ε W (y) and b are upper- bounds of E y , so is ε W (y) ∧ b. Hence, ε W (y) ∧ b = E y = α W (y)
. By a similar reasoning as in Section 2.2, we get the following result.

Proposition 2 Let W ∈ M n and ε W be the function defined by [START_REF] Blusseau | Tropical and morphological operators for signal processing on graphs[END_REF]. Then ε W is an erosion mapping L to L if and only if W is column-0-astic.

If W is also row-0-astic, then ε W = α W is the adjoint of δ W , as stated next.

Proposition 3 Let W ∈ M n and δ W and ε W be the functions defined by ( 1) and ( 2), respectively. Then (ε W , δ W ) is an adjunction on L if and only if W is doubly-0-astic. Furthermore, (ε W , δ W ) is an adjunction on L with δ W extensive (and ε W anti-extensive) if and only if W is a CMW matrix.

Proof. Most of the points have already been addressed above or are straightforward from Proposition 1. To see that δ W extensive implies w ii = 0 for all i, just remark that w ii < 0 would imply δ W (e (i)

) i < b = e (i) i . ⊓ ⊔ 2.4 Generality of (ε W , δ W )
The dilation δ W , already introduced in [START_REF] Blusseau | Tropical and morphological operators for signal processing on graphs[END_REF][START_REF] Maragos | Chapter Two -Representations for morphological image operators and analogies with linear operators[END_REF][START_REF] Velasco-Forero | Nonlinear operators on graphs via stacks[END_REF], can be viewed as a generalisation of the non-local and adaptive mathematical morphology [START_REF] Salembier | Study on nonlocal morphological operators[END_REF][START_REF] Velasco-Forero | On nonlocal mathematical morphology[END_REF] on signals and images. Each column W :,j of W represents the structuring function corresponding to pixel (or instant) j.

As pointed out in [START_REF] Heijmans | The algebraic basis of mathematical morphology I. dilations and erosions[END_REF][START_REF] Maragos | Chapter Two -Representations for morphological image operators and analogies with linear operators[END_REF], the dilations that can be written as matrix-based max-plus products like Eq. ( 1) are the shift (or vertical-translation) invariant ones. However the result stated in [START_REF] Heijmans | The algebraic basis of mathematical morphology I. dilations and erosions[END_REF][START_REF] Maragos | Chapter Two -Representations for morphological image operators and analogies with linear operators[END_REF] does not directly apply to our setting where the lattice L is different from the lattice of scalars which define vertical translation of signal values, usually R ∪ {-∞, +∞}. Still, the same idea holds here with some adaptation, as stated in the next proposition.

Proposition 4 Let δ : L → L be a dilation. Then there exists W ∈ M n such that δ = δ W if and only if ∀λ ≤ 0, ∀x ∈ L, δ (λ + x) ∨ a = λ + δ(x) ∨ a. (3)
In that case, the matrix W whose j-th column is W :,j = δ(e (j) )-b for 1 ≤ j ≤ n, is such a representing matrix.

We see that this class of dilations is very broad and covers the most commonly used in morphological image and signal processing: dilations based on structuring functions, possibly non-local, varying in space and non-flat.

Proof (Proposition 4). If δ = δ W for some W ∈ M n , then it is straightforward to check that δ verifies Eq. ( 3). Conversely, suppose δ verifies Eq. ( 3). Then we first remark that δ(b) = b. Indeed, on the one hand, δ(a) = a as a = L and δ is a dilation mapping L to L. On the other hand, δ(a

) = δ (a -b) + b = (a -b) + δ(b) by Eq. (3). Hence a = (a -b) + δ(b) which means that δ(b) = b.
As a consequence: for any i ∈ 1, n , there is a j i ∈ 1, n such that δ(e (j i ) ) i = b. This is simply because b = 1≤j≤n e (j) so b = δ 1≤j≤n e (j) = 1≤j≤n δ(e (j) ), which means that, for any i, b = 1≤j≤n δ(e (j) ) i and finally that δ(e (j i ) ) i = b for some j i , as the supremum is reached here. Now, let x ∈ L. Then it can be decomposed as x = 1≤j≤n (λ j + e (j) ) ∨ a with λ j = x jb ≤ 0. Hence, as δ is a dilation verifying Eq. ( 3), we get δ(x) = 1≤j≤n λ j + δ(e (j) ) ∨ a . We now use the result stated just above: for any i ∈ 1, n there is a

j i ∈ 1, n such that λ ji + δ(e (j i ) ) i = x ji -b + b = x ji ≥ a.
Therefore, 1≤j≤n λ j + δ(e (j) ) ∨ a = 1≤j≤n λ j + δ(e (j) ) from which we finally get

δ(x) = 1≤j≤n λ j + δ(e (j) ) = 1≤j≤n (x j -b) + δ(e (j) ) = 1≤j≤n x j + [δ(e (j) ) -b] (4)
which is exactly W ⊗ x for W the matrix with columns W :,j = δ(e (j) )b for 1 ≤ j ≤ n.

⊓ ⊔

Note that the dual of Proposition 4 obviously holds: the erosions ε : L → L which can be written as ε W for some W ∈ M n are those for which

∀λ ≥ 0, ∀x ∈ L, ε (λ + x) ∧ b = λ + ε(x) ∧ b. (5)
To show this it is sufficient to see that ε verifies (5) if and only if the dilation δ = ε(• c ) c verifies (3), and recall that δ

W (• c ) c = ε W T (•).

Equivalent dilations and erosions

In Proposition 4 we exhibited one possible matrix W ∈ M n that represents a dilation, but this matrix is not unique. In this section we characterise the set of such matrices and show that it is a complete lattice. Since we are interested in adjunctions (ε W , δ W ), following Proposition 3 we focus on the set of matrices in M n that are doubly-0-astic, which we denote by D 0 (n). Let the equivalence relation defined for any two matrices

A, B ∈ D 0 (n) by A ∼ B ⇐⇒ δ A = δ B ⇐⇒ ∀x ∈ L, δ A (x) = δ B (x) (6)
and note

C W = {M ∈ D 0 (n), M ∼ W } the equivalence class of any W ∈ D 0 (n).
We provide an easy characterisation of C W that will show useful in numerical computations of the morphological operators defined earlier. For any u ∈ R max let I u denote the matrix in M n whose coefficients are all equal to u. Then two equivalent matrices are characterised as follows.

Proposition 5 Let M, W ∈ D 0 (n). Then M ∈ C W ⇐⇒ M ∨ I a-b = W ∨ I a-b ⇐⇒ m ij = w ij if w ij > a -b m ij ≤ a -b otherwise. (7)
This means that if W has coefficients not larger than ab, these can be set to any value not larger than ab, including -∞, and can therefore be ignored in the computation of δ W (x).

Proof (Proposition 5). The second equivalence is just a matter of writing, so we prove the first one. Let us first notice that for any x ∈ L,

I (a-b) ⊗ x ≤ a. Therefore ∀x ∈ L, (W ∨ I (a-b) ) ⊗ x = (W ⊗ x) ∨ (I (a-b) ⊗ x) = W ⊗ x, since W ⊗ x ≥ a, and this holds for M too. Hence, if M ∨ I a-b = W ∨ I a-b , then for any x ∈ L, W ⊗ x = (W ∨ I a-b ) ⊗ x = (M ∨ I a-b ) ⊗ x = M ⊗ x, which means M ∈ C W .
Conversely, suppose that M ∼ W and that w i0j0 > ab for some i 0 , j 0 ∈ 1, n . Let x = e (j0) ∈ L, i.e. x j0 = b and x j = a ∀j = j 0 . The 0-asticity of W and

M implies (W ⊗ x) i0 = b + w i0j0 and (M ⊗ x) i0 = b + m i0j0 , hence m i0j0 = w i0j0 . We have just shown that ∀i, j ∈ 1, n , (w ij > a -b ⇒ w ij = m ij )
and by symmetry of the equivalence relation (

m ij > a -b ⇒ w ij = m ij ), which combined yields max(m ij , a -b) = max(w ij , a -b). So finally M ∼ W ⇒ M ∨ I a-b = W ∨ I a-b . ⊓ ⊔ While it is clear that if A, B ∈ C W then A ∨ B ∈ C W , the characterisation in Proposition 5
shows that C W is also closed under infimum, that is:

A ∧ B ∈ C W .
This has the following straightforward consequence.

Proposition 6 Let W ∈ D 0 (n) and ≤ the partial ordering on C W defined by

A ≤ B ⇐⇒ A ∨ B = B ⇐⇒ a ij ≤ b ij ∀i, j ∈ 1, n . Then -(C W , ≤
) is a complete lattice (with coefficient-wise supremum and infimum);

-Its greatest element is W = W ∨ I a-b ; -Its smallest element is W , defined by w ij = w ij if w ij > a -b -∞ otherwise.

Iterated operators and granulometries

In this section, given W ∈ D 0 (n) and p ∈ N * , we focus on the iterated dilations and erosions δ p W and ε p W , as well as their sup and inf integrations, that we note respectively D W ) are adjunctions. We note respectively γ

[p] W =δ p W ε p W and G [p] W =D [p] W E [p]
W their corresponding openings.

Note that if δ W is extensive, or equivalently if W is a CMW matrix (Prop. 3), then these adjunctions are equal:

(ε p W , δ p W ) = (E [p] W , D [p]
W ). As this is not true in general, both adjunctions are worth studying. In particular, we shall examine whether (γ

[p] W ) p∈N * and (G [p]
W ) p∈N * define granulometries, that is to say families of openings that are decreasing with p. The answer is yes and it is a general result that does not depend on the representation of the adjunction.

Proposition 7 Let (ε, δ) be an adjunction on a complete lattice. For any integer p ∈ N * , let us note γ p = δ p ε p and G p = D p E p the openings associated to the adjunctions (ε p , δ p ) and E p = 1≤k≤p ε k , D p = 1≤k≤p δ k , respectively. Then (γ p ) p∈N * and (G p ) p∈N * are granulometries.

Proof. We first show that the family of openings (γ p ) p≥1 decreases with p, hence a granulometry. This is straightforward by writing γ p+1 = δ p+1 ε p+1 = δ p γ 1 ε p ≤ δ p ε p = γ p . Secondly, regarding (G p ) p≥1 , we show G p+1 ≤ G p by proving that G p G p+1 = G p+1 . We obtain this by remarking that D p+1 = D p (id δ), which makes it an invariant of G p :

G p D p+1 = D p E p D p (id δ) = D p (id δ) = D p+1 . Then we can conclude G p G p+1 = G p D p+1 E p+1 = D p+1 E p+1 = G p+1 . ⊓ ⊔
To conclude this section, let us write δ p W , ε p W , D

W and E

[p]

W as dilations and erosions represented by one suitable doubly-0-astic matrix. This will help in their graph interpretation of the next section. The associativity of ⊗ yields ∀x ∈ L,

δ p W (x) = W ⊗ . . . ⊗ W ⊗ x = W p ⊗ x, therefore δ p W = δ W p . We obtain similarly ε p W = ε W p . The distributivity of ⊗ over ∨ yields D [p] W (x) = p k=1 δ k W (x) = p k=1 (W k ⊗ x) = ( p k=1 W k ) ⊗ x therefore D [p] W = δ Sp(W ) , with S p (W ) = p k=1 W k . Similarly, E [p] W = ε Sp(W )
. Note that by the same arguments and Proposition 3, we get that D 0 (n) is closed under ⊗ and ∨.

Graph interpretations

Weighted graphs

Let W ∈ M n and G(W ) = (V, E) be a weighted and directed graph containing n vertices whose n × n adjacency matrix is W , with the convention that w ij > -∞ if and only if (i, j) ∈ E. We now recall that a path from vertex i to vertex j in G(W ) is a tuple π = (k 1 , . . . , k l ) of vertices such that k 1 = i, k l = j, and (k p , k p+1 ) ∈ E for 1 ≤ p ≤ l -1. The length of the path, denoted by ℓ(π), is l -1 (the number of its edges). For p ≥ 1, Γ (p) ij (W ) denotes the set of paths from i to j in G(W ) of length p and Γ (∞) ij (W ) the set of all paths from i to j. The weight of a path π = (k 1 , . . . , k l ), denoted by ω(π), is the sum ω(π) = l-1 p=1 w kpkp+1 .

Iterated operators

Recall that for W ∈ M n and p ∈ N * , W p is the p-th power of W in the ⊗ sense, and S p (W ) is the matrix defined in Section 2.6, denoted by S p here for simplicity. We note respectively w (p) ij and s

[p] ij their coefficients. The following result is well known in tropical algebra and graph theory [START_REF] Carré | An algebra for network routing problems[END_REF][START_REF] Cuninghame-Green | Minimax algebra[END_REF], and will help interpret the operators defined earlier. It can be proved by induction.

Proposition 8 Let W ∈ M n and p ∈ N * . Then for any 1 ≤ i, j ≤ n, w

(p) ij = max ω(π), π ∈ Γ (p) ij (W ) and s [p] ij = max    ω(π), π ∈ 1≤k≤p Γ (k) ij (W )    (8) with the convention max(∅) = -∞.
The equations in ( 8) are equivalent to saying that ij ) is the maximal weight over the set of paths from vertex i to vertex j of length exactly (resp. at most) p.

1. w (p) ij > -∞ (resp. s [p] ij > -∞) if
Therefore the graphs G(W p ) and G(S p ) have the same set of vertices as the original graph G(W ), but an edge exists between vertices i and j in G(W p ) (resp. G(S p )) whenever there is a path of length exactly (resp. at most) p from i to j in G(W ). The weights associated with this new edge are the maximal weights over the corresponding set of paths. Now if W ∈ D 0 (n), following Section 2.6 we get, for x ∈ L and i ∈ 1, n :

δ p W (x) i = j∈N p i {x j + w (p) ij } , ε p W (x) i = j∈ Ň p i {x j -w (p) ji } (9)
and D

[p]

W (x) i = j∈N p i {x j + s [p] ij } , E [p] W (x) i = j∈ Ň k i {x j -s [p] ji } (10)
where N p i is the set of neighbours of vertex i in G(W p ) or, equivalently, the set of vertices in G(W ) that can be reached from i through a path of length

p; Ň p i = j ∈ {1, . . . , n}, i ∈ N p j ; N p i = ∪ 1≤k≤p N k i and Ň p i = ∪ 1≤k≤p Ň k i .
Hence these dilations and erosions are suprema and infima of "penalised" values over extended neighbourhoods induced by the original graph. The penalization is given by the strength of the connection between vertices: the closer the penalising weight to zero, the more the neighbours' value contributes to the result. The fact that we can restrict the supremum and infimum over graph neighbourhoods in ( 9) and ( 10) is due to the weight values being -∞ outside these neighbourhoods, hence not contributing to the supremum and infimum.

Path interpretation of the opening G [p] W

The goal of this section is to show that G

[p]

W can be interpreted similarly to a path opening [START_REF] Heijmans | Path openings and closings[END_REF], in the sense that it preserves bright values that are connected to other bright values forming long enough paths in a graph. We can first remark that for any x ∈ L, i ∈ 1, n and t ∈ [a, b]:

G [p] W (x) i ≥ t ⇐⇒ ∃j ∈ N p i , such that ∀l ∈ Ň p j x l ≥ t -s [p] ij + s [p] lj , (11)
which is straightforward from the expressions in [START_REF] Lerallut | Image filtering using morphological amoebas[END_REF], as

G [p] W = D [p] W E [p]
W . This directly yields

G [p] W (x) i = t ∈ [a, b], ∃j ∈ N p i , ∀l ∈ Ň p j , x l ≥ t -s [p] ij + s [p] lj . (12)
In the case of binary weights, i.e. w ij = 0 if vertex j is neighbour of i in G and w ij = -∞ otherwise, which corresponds to a non-weighted graph, then s [START_REF] Maragos | Chapter Two -Representations for morphological image operators and analogies with linear operators[END_REF] and [START_REF] Ronse | Why mathematical morphology needs complete lattices[END_REF]. Therefore, if G

[p] ij = s [p] lj = 0 in
[p] W (x) i ≥ t, then there is a vertex j which is at most p steps away from i, such that all paths of length at most p and ending in j, including those of length exactly p and passing through i (if they exist), show values larger than t. In the general case, the additional term -s

[p] ij +s [p]
lj modulates this constraint in function of the strength of the connection of i and the other vertices of Ňj , to j.

Links to the max-plus spectral theory

Now we present the consequences and interpretations of some results from the spectral theory in max-plus algebra. We first report definitions from [START_REF] Cuninghame-Green | Minimax algebra[END_REF] necessary to Theorem 1 (also from [START_REF] Cuninghame-Green | Minimax algebra[END_REF]). Then we draw the links to our setting and more particularly in the case of a symmetric matrix, corresponding to a non-directed graph. In all this section, W ∈ M n .

General definitions and results

Definition 3 (Eigenvector, eigenvalue [START_REF] Cuninghame-Green | Minimax algebra[END_REF]) Let x ∈ R n max and λ ∈ R max . Then x is an eigenvector of W with λ as corresponding eigenvalue if W ⊗ x = λ ⊗ x = λ + x. If there exists finite x and λ solutions to this equation, we say that the eigenproblem is finitely soluble.

In the graph G(W ), a path (k 1 , . . . , k l ) is called a circuit if k 1 = k l . We will note C(W ) the set of all circuits of G(W ). Circuits allow us to distinguish another class of matrices in M n , called definite They are important to the present framework as they include the doubly-0-astic matrices.

Definition 4 (Definite matrix [START_REF] Cuninghame-Green | Minimax algebra[END_REF]) W is said definite if max c∈C(W ) ω(c) = 0. In other words, all the circuits of G(W ) have non positive weights, and at least one circuit c * , called a zero-weight circuit, achieves ω(c * ) = 0.

To see that if W is row or column-0-astic, then it is definite, it is sufficient to build an increasing path with zero-weight, until one vertex repeats. The path can be initialized with any vertex j 1 . Then given the current path (j 1 , . . . , j m ), we extend it by adding a vertex j m+1 such that w jmjm+1 = 0. This is always possible thanks to the row or column-0-asticity of W . Since there are n distinct vertices in G(W ), an index will repeat after at most n iterations.

Definition 5 (Eigen-node, equivalent eigen-nodes [START_REF] Cuninghame-Green | Minimax algebra[END_REF]) Let W be a definite matrix. An eigen-node is any vertex in G(W ) belonging to a zero-weight circuit. Two eigen-nodes are said equivalent if there is a zero-weight circuit passing through both of them.

In [START_REF] Cuninghame-Green | Minimax algebra[END_REF], S n (W ) = 1≤k≤n W k is denoted by ∆(W ) and called the metric matrix. Recall that for i, j ∈ 1, n , ∆(W ) ij is the maximal weight over the set of paths from vertex i to vertex j of length at most n, in G(W ) (Prop. 8). If W is definite, circuits have non-positive weights in G(W ) and therefore any path longer than n can be reduced to a shorter path with non larger weight. Hence, ∆(W ) ij is actually the maximal weight over the set of all paths from i to j. This provides an easy characterisation of eigen-nodes for W definite: j is an eigen-node of G(W ) if and only if ∆(W ) jj = 0. Furthermore, the j-th column ξ j of ∆(W ) is a map of the ancestors of j in G(W ). It tells which vertices can reach j and at which cost. Definition 6 (Fundamental eigenvectors, eigenspace [START_REF] Cuninghame-Green | Minimax algebra[END_REF]) Let W be a definite matrix. Then a fundamental eigenvector of W is any j-th column ξ j of ∆(W ), where j is an eigen-node. Two fundamental eigenvectors are said equivalent if their associated eigen-nodes are equivalent (see Definition 5). Let E = {ξ i1 , ξ i2 , . . . , ξ i k } be a set of k ≥ 1 fundamental eigenvectors of W , all pairwise non-equivalent. The set E is said to be a maximal set of non-equivalent fundamental eigenvectors if any other fundamental eigenvector of W is equivalent to one of the eigenvectors in E. In this case the set { k j=1 x j + ξ ij , x ∈ R k max } is called the eigenspace of W and does not depend on E (see [START_REF] Cuninghame-Green | Minimax algebra[END_REF], Lemma 24-1).

Theorem 1 ([6]

) Let W be a doubly-0-astic (hence definite) matrix. Then the following statements are valid:

-For any fundamental eigenvector ξ j of W (finite or not), W ⊗ ξ j = ξ j .

-The eigenproblem is finitely soluble.

-If two fundamental eigenvectors are equivalent, then they are equal.

-Any finite eigenvector is associated to the eigenvalue λ = 0, and lies in the eigenspace of W .

Consequences and interpretations

In general. As said, the results of the previous section apply to our setting since we consider adjunctions represented by doubly-0-astic matricesa. For W ∈ D 0 (n), ∆(W ) is also in D 0 (n) and the corresponding opening δ

∆(W ) ε ∆(W ) is G [n]
W . By definition, G

[n] W (x) projects x ∈ L onto δ ∆(W ) (L), which is the set { n j=1 y j + ξ j , y ∈ L} of max-plus combinations of columns of ∆(W ). Theorem 1 tells that this decomposition can be split as G

[n] W (x) = u ∨ v, where u lies in the eigenspace of W and v is a max-plus combination of the ξ j s which are not fundamental eigenvectors. This decomposition may be sparser than the original one, as the dimension of the eigenspace of W , i.e. Card(E), can be lower than the number of fundamental eigenvectors.

The case of symmetric W ∈ D 0 (n). This case corresponds to considering a nondirected graph supporting the signal x. As the adjacency relationship is often based on a symmetrical function on pairs of vertex values, this assumption covers many practical cases (e.g. [START_REF] Blusseau | Tropical and morphological operators for signal processing on graphs[END_REF][START_REF] Blusseau | Adaptive Anisotropic Morphological Filtering Based on Co-Circularity of Local Orientations[END_REF]). The main consequence of W ∈ D 0 (n) symmetric is that every vertex j is an eigen-node: for any j ∈ 1, n there is i such that w ij = 0 = w ji and therefore (j, i, j) is a zero-weight circuit. This entails three other consequences.

First, ∆(W ) jj = 0 for every j ∈ 1, n , following the characterisation of eigen-nodes described earlier, which implies that δ ∆(W ) = D W anti-extensive (Prop. 3). Secondly, W ⊗ ξ j = ξ j for every column ξ j of ∆(W ), which implies W k ⊗ ξ j = ξ j for 1 ≤ k ≤ n, hence ∆(W ) ⊗ ξ j = ξ j and finally ∆(W ) ⊗ ∆(W ) = ∆(W ). This means D W is exactly the eigenspace of W .

When W is symmetric, a maximal set of k non-equivalent fundamental eigenvectors {ξ i1 , ξ i2 , . . . , ξ i k }, k ≤ n, can be seen as negative distance maps to the k corresponding eigen-nodes G(W ), as they contain the optimal cost (maximal weight) between any vertex and the eigen-nodes4 . Hence we can picture the aspect of G

[n] W (x), for x ∈ L: it is the upper-envelope of the largest vertical translations of these distance maps that are dominated by x. Therefore, adapting G(W ) to x by well connecting vertices within relevant structures preserves these structures under the filter G

[n] W , as shown in [START_REF] Blusseau | Tropical and morphological operators for signal processing on graphs[END_REF][START_REF] Blusseau | Adaptive Anisotropic Morphological Filtering Based on Co-Circularity of Local Orientations[END_REF]. In practice, n might be large, such as the number of pixels of an image. Since (G

[p]

W ) 1≤p≤n is a granulometry, we know that G

[n]

W can be approximated by G

[p]

W with increasing p.

Conclusion

In this paper we consolidated the basis of the representation of adjunctions by matrices in max-plus algebra. We showed that it is a very flexible framework that generalises many types of morphological adjunctions. In particular, it allows describing precisely the behaviour of iterated operators based on spatially-variant, non-flat structuring functions. This is made possible by their graph interpretation and spectral results in max-plus algebra. Future works shall investigate further the insights that max-plus algebra can bring to mathematical morphology through this framework.

W

  is extensive andε ∆(W ) = E [n]

W

  are idempotent. They are therefore a closing and an opening respectively and E an adjunction (ε, δ) for which ε is an opening and δ a closing verifies ε = δε (and δ = εδ). The third consequence is the following.Corollary 1 If W ∈ D 0 (n) is symmetric, then the set of invariants of G [n]

  and only if there is at least a path in G(W ) from vertex i to vertex j of length exactly (resp. at most) p; 2. w

	(p) ij (resp. s	[p]

An online demo for[START_REF] Blusseau | Adaptive Anisotropic Morphological Filtering Based on Co-Circularity of Local Orientations[END_REF] is available: https://bit.ly/anisop demo.

Note that ∆(W ) is a metric, not exactly between vertices, but between their equivalence classes induced by Def

. 5, as all vertices are eigen-nodes when W is symmetric.


S Mazzi

J Garcia

D Zarzoso

Ye O Kazakov

J Ongena

M Dreval

M Nocente

Ž Štancar

G Szepesi

J Eriksson

A Sahlberg

S Benkadda

Enhanced performance in fusion plasmas through turbulence suppression by MeV ions

Alpha particles with energies of the order of MeV will be the main source of plasma heating in future magnetic confinement fusion reactors. Instead of heating fuel ions, most of the energy of alpha particles is transferred to electrons in the plasma. Furthermore, alpha particles can also excite Alfvénic instabilities, which were previously considered as detrimental for the performance of the fusion device. Here we report improved thermal ion confinement in the presence of MeV ions and strong fast ion driven Alfvénic instabilities in recent experiments on the Joint European Torus. Detailed transport analysis of these experiments reveals turbulence suppression through a complex multi-scale mechanism that generates large-scale zonal flows. This holds promise for a more economical operation of fusion reactors with dominant alpha particle heating and ultimately cheaper fusion electricity.

The urgent need for environmentally friendly energy sources becomes increasingly more important. Among several options, the fusion reaction between the hydrogen isotopes deuterium (D) and tritium (T) holds the promise for a safe clean and almost inexhaustible energy production: D + T → 4 He (3.5 MeV) + n (14.1 MeV), with an alpha particle ( 4 He) and a neutron (n) as fusion products. This reaction needs temperatures around 150 million K and can be realized in magnetic-confinement fusion devices [1]. The largest fusion device currently in operation is the Joint European Torus (JET) [2]. The successor device ITER aims to demonstrate the scientific and technological matureness of magnetic confinement fusion [3].

Both JET and ITER are based on the tokamak concept, a toroidal configuration schematically shown in Figure 1. Large temperature gradients (∼100 million K per meter) are unavoidably present in these devices because of the required high plasma temperature in the center and the necessity for a cold plasma edge. These gradients create instabilities and turbulence, often with very different spatio-temporal scale lengths [4]. One of the most important is the ion temperature gradient (ITG) instability with a characteristic scale length on the order of the ion gyroradius, usually much smaller than the plasma size (e.g., ∼10 -3 m vs. ∼1 m in JET). The microscopic ITG instability largely limits plasma temperatures that can be achieved in a fusion device [5].

The success of magnetic confinement fusion as an energy source relies crucially on reaching high temperatures for the D and T ions. Fusion-born MeV-range alpha particles are the main source of central plasma heating in ITER and future fusion power plants. Yet, these highly energetic ions heat primarily electrons rather than thermal ions through Coulomb collisions. The extrapolation of plasma heating by alpha particles in future devices is not straightforward partly because of the impact of fast ions on plasma turbulence. Recent progress in this developing field of research showed that ion scale turbulence can be partly reduced in the presence of fast ions with energies of ∼100 keV [6][7][8][START_REF] Di Siena | Electromagnetic turbulence suppression by energetic particle driven modes[END_REF]. However, at these energies fast ions provide dominant bulk ion heating in strong contrast to electron heating from alphas in a fusion reactor.

Alpha particles can also excite a range of instabilities, which are often considered detrimental for plasma confinement. Modes driven by fast ions can cause increased radial transport of energetic ions [START_REF] Fasoli | Chapter 5: Physics of energetic ions[END_REF][START_REF] Gorelenkov | Energetic particle physics in fusion research in preparation for burning plasma experiments[END_REF][START_REF] Todo | Introduction to the interaction between energetic particles and Alfvén eigenmodes in toroidal plasmas[END_REF][START_REF] Heidbrink | Mechanisms of energetic-particle transport in magnetically confined plasmas[END_REF] and enhance plasma turbulence [START_REF] Zarzoso | Impact of Energetic-Particle-Driven Geodesic Acoustic Modes on Turbulence[END_REF]. Among possible instabilities driven by alpha particles in ITER, toroidicity-induced Alfvén eigenmodes (TAEs) [START_REF] Cheng | Low-n shear Alfvén spectra in axisymmetric toroidal plasmas[END_REF][START_REF] Fu | Excitation of the toroidicity-induced shear Alfvén eigenmode by fusion alpha particles in an ignited tokamak[END_REF] are of particular concern [START_REF] Pinches | Energetic ions in ITER plasmas[END_REF]. Yet numerical analysis in [START_REF] Di Siena | Electromagnetic turbulence suppression by energetic particle driven modes[END_REF] showed that in plasmas with ion heating from ∼100 keV ions, i.e. under conditions far from those in a fusion reactor, reduced turbulence can be obtained when TAEs are marginally stable, leaving the question open on whether improved thermal ion confinement can be reached experimentally in plasmas with MeV-range fast ions and fully destabilized TAEs.

Thermal ion energy fluxes perpendicular to the magnetic surfaces can be reduced because of the appearance of intense poloidally directed shear flows, known as zonal flows [START_REF] Diamond | Zonal flows in plasmas -a review[END_REF]. In fact, zonal flows, which are solely generated by nonlinear interactions, are a common phenomenon in nature and can be remarkably stable as, e.g., the famous belts of Jupiter [START_REF] Heimpel | Simulation of equatorial and high-latitude jets on Jupiter in a deep convection model[END_REF]. The generation of zonal flows by energetic ions and fast ion driven modes, in particular TAEs, was put forward in simplified contexts in several theoretical studies recently [START_REF] Todo | Nonlinear magnetohydrodynamic effects on Alfvén eigenmode evolution and zonal flow generation[END_REF][START_REF] Qiu | Effects of energetic particles on zonal flow generation by toroidal Alfvén eigenmode[END_REF]. As the dynamics of ITER plasmas will be determined in part by the nonlinear interaction between MeV-range alphas, Alfvén eigenmodes destabilized by these alphas, and microturbulence, studying this complex interplay at ITER-relevant conditions becomes an urgent necessity for future progress in fusion science.

In this paper, we identify reproducible conditions and provide experimental evidence for the suppression of ion-scale turbulence in fusion-reactor-grade plasmas with strong electron heating from MeV-range fast ions and simultaneously destabilized Alfvén eigenmodes. Corroborated by detailed turbulence analysis, our results pave the way towards enhanced performance of future fusion reactors with strong alpha particle heating.

Improved thermal ion confinement in plasmas with MeV ions

The impact of fast ions on plasma dynamics was studied in JET D- 3 He plasmas (n( 3 He)/n e ≈ 20-30%). While the presence of TAEs is usually accompanied by a loss of energy and particle confinement [START_REF] Fasoli | Chapter 5: Physics of energetic ions[END_REF][START_REF] Gorelenkov | Energetic particle physics in fusion research in preparation for burning plasma experiments[END_REF][START_REF] Chen | Physics of Alfvén waves and energetic particles in burning plasmas[END_REF][START_REF] Breizman | Major minority: energetic particles in fusion plasmas[END_REF][START_REF] Heidbrink | Basic physics of Alfvén instabilities driven by energetic particles in toroidally confined plasmas[END_REF], here we show that reduced ion heat losses and high ion temperatures can be reached in plasmas with MeV-range fast ions and TAEs driven by fast ions. Figure 2 compares JET pulses #94704 (∼100 keV ions; blue lines) and #94701 (MeV-range ions; red lines) that were performed at the same operational parameters (L-mode, B 0 = 3.7 T, I p = 2.5 MA, n e0 ≈ 6 × 10 19 m -3 ) and total auxiliary heating power P aux = 14 MW, but differing in the characteristics of the fast ion population.

In pulse #94704 (blue lines in Fig. 2), Neutral Beam Injection (NBI) was the only heating system, providing fast D ions with energies up to 100 keV. These moderately energetic ions deposit most of their energy to bulk ions, resulting in plasmas with T i /T e > 1 (e.g., T i /T e ≈ 1.4 measured at ρ tor = 0.2). The corresponding electron and ion temperature profiles at P aux = 14 MW (t = 9.0 s) are shown in Figs. 2(f) and (g). The injected fast NBI ions are sub-Alfvénic (v NBI /v A ≈ 0.4) and no Alfvén eigenmodes were observed in this plasma.

In the comparison pulse #94701 (red lines in Fig. 2), fast D ions from NBI were accelerated to much higher energies (up to ∼2-3 MeV; see Fig. S1 in the Supplementary Information) with waves in the ion cyclotron range of frequencies (ICRF) using the three-ion D-(D NBI )-3 He scheme (P NBI = 8 MW, P ICRF = 6 MW) [START_REF] Kazakov | Efficient generation of energetic ions in multi-ion plasmas by radio-frequency heating[END_REF][START_REF] Ongena | Synergetic heating of D-NBI ions in the vicinity of the mode conversion layer in H-D plasmas in JET with the ITER like wall[END_REF][START_REF] Nocente | Generation and observation of fast deuterium ions and fusion-born alpha particles in JET D-3 He plasmas with the 3-ion radio-frequency heating scenario[END_REF][START_REF] Kazakov | Physics and applications of three-ion ICRF scenarios for fusion research[END_REF]. The presence of MeV-range ions resulted simultaneously in long-period sawteeth, destabilization of various Alfvén eigenmodes, as well as in the increased plasma stored energy and D-D neutron rate. As a result of the collisional slowing-down, these highly energetic deuterons heat predominantly electrons (similar to alpha particle heating in future fusion reactors), giving rise to a strongly peaked electron temperature profile with much larger T e0 , see Fig. 2(f). Interestingly, the ion temperature in this core-electron-heated plasma was at least as high as in the NBI-only pulse #94704 with more ion heating and T i /T e ≈ 1 was reached.

We focus on the transport characteristics in the central region of these plasmas (ρ tor 0.4), where most of the alpha particles will be generated in a future fusion reactor. The electron and thermal ion heat diffusivities computed with TRANSP [START_REF] Hawryluk | An empirical approach to tokamak transport[END_REF][START_REF] Ongena | Numerical Transport Codes[END_REF], see Figures 2(h) and (i), show that ions are the primary heat loss channel in the NBI-only pulse #94704 (χ i > χ e ): e.g., at ρ tor = 0.23, χ i ≈ 1.9 m 2 s -1 and χ e ≈ 0.6 m 2 s -1 . In contrast, in pulse #94701 with highly energetic ions, χ i ≈ 0.8 m 2 s -1 is much reduced at this radial location, while χ e remains almost unchanged. The comparison of two χ i -profiles shown in Fig. 2(i) clearly indicates a significant improvement of thermal ion confinement in the plasma core of the pulse with electron heating by fast ions. Note that in the absence of ion temperature measurements at ρ tor < 0.2, the near-axis χ i were computed using a global third-order polynomial fit for T i (routinely employed in TRANSP modelling of JET experiments), or assuming T i = T e ; see the Methods section.

As a next step, predictive modeling of T i in pulse #94701 was undertaken with the CRONOS code [START_REF] Artaud | The CRONOS suite of codes for integrated tokamak modelling[END_REF], using all relevant parameters, including measured T e and n e , and heat source profiles calculated by TRANSP as input. Two different χ i -profiles were assumed in this analysis, see Fig. 2(i). The impact of the reduced ion heat losses on the achieved T i in #94701 is clearly illustrated by comparing the two predicted T i -profiles shown in Fig. 2(g). The red dotted line illustrates the predicted ion temperature in #94701 using the TRANSP-computed χ i for this pulse. The good agreement of the predicted and measured T i profiles verifies the consistency of CRONOS and this synthetic modelling approach. This analysis also shows that a much lower T i would have been reached in #94701 (red dashed line), if the thermal ion heat diffusivity were the same as in the NBI-only pulse #94704. However, the experimentally observed T i in #94701 is significantly larger than the latter T i prediction, another confirmation for reduced ion heat losses in the presence of MeV ions.

The well-known mechanisms stabilizing ITG turbulence such as the ratio T i /T e and the E × B shearing rate (γ E×B ) cannot explain the lower χ i in the core region ρ tor < 0.4 for pulse #94701. Indeed, both parameters are lower in this pulse as compared to those in the NBI-only pulse #94704: e.g., at ρ tor = 0.23, T i /T e ≈ 1 vs. 1.4 and

γ E×B ≈ 1.4 m s -1 vs. 2.2 m s -1 .
The two comparison pulses are not only inherently different in the heating channel for fast ions (dominant ion vs. electron heating), but also feature an essential difference in the MHD behavior. The sub-Alfvénic NBI ions in #94704 did not destabilize any Alfvén eigenmodes. In contrast, a variety of TAE modes in the frequency range ∼200 kHz with toroidal mode numbers n = 2 -6 were driven by high-energy ions in pulse #94701, see Fig. 2(e). The observed mode frequencies are in good agreement with the theoretically predicted frequencies of TAEs, f

(lab) TAE = f (plasma) TAE + nf rot . In this formula, f (plasma) TAE = v A0 /(4πq TAE R 0 ) ≈ 166 -191 kHz, q TAE 1 + 1/(2n) ≈ 1.1 -1.3 [24],
v A0 is the Alfvén velocity in the plasma center, and f rot ≈ 7 -8 kHz is the local plasma rotation frequency, as measured by the CXRS system. The localization of TAEs outside the q = 1 surface in this ICRF+NBI plasmas has been independently confirmed by correlation reflectometer measurements, indicating that R TAE ≈ 3.22 -3.36 m (see Fig. S3(a) in the Supplementary Information).

The correlation reflectometer also provides local density fluctuation spectra, used to characterize the turbulence characteristics in the plasma. The performed analysis at R = 3.25 m and R = 3.3 m (see Figs. S3(b) and (c) in the Supplementary Information) highlights that the lowest fluctuation amplitudes were reached in plasmas with MeV-range fast ions and destabilized TAEs. Furthermore, a clear reduction of the thermal ion heat diffusivity was observed at these radial locations, see Fig. 2(i). As the characteristics of the fast ions in pulse #94701 (T fast /T e ≈ 34 and n fast /n e ≈ 3%) are close to the predicted conditions for alpha particles in ITER [START_REF] Garcia | Isotope and fast ions turbulence suppression effects: Consequences for high-β ITER plasmas[END_REF], the identified mechanism holds promise for optimizing the performance of future fusion reactors with strong alpha-particle heating.

Note that the observed sawtooth dynamics in pulse #94701 is similar to earlier JET experiments with ICRF and monster sawteeth reported in [START_REF] Campbell | Stabilization of sawteeth with additional heating in the JET tokamak[END_REF][START_REF] Graves | Experimental verification of sawtooth control by energetic particles in ion cyclotron resonance heated JET tokamak plasmas[END_REF]. However, in strong contrast to [START_REF] Campbell | Stabilization of sawteeth with additional heating in the JET tokamak[END_REF], a much smaller plasma volume was enclosed within the sawtooth inversion radius in the JET experiments reported here, thereby limiting the impact of the sawtooth stabilization itself on the improved plasma confinement. We also note that the power balance analysis used in this paper to infer the thermal conductivities ignores any possible anomalous ion heating effects such as, e.g., reported in [START_REF] Gates | Ion heating by fast-particle-induced Alfvén turbulence[END_REF]. In the absence of reliable measurements of high-frequency Alfvén eigenmodes (> 1 MHz), their potential impact in the reported plasmas is difficult to assess.

Numerical modeling of turbulence suppression

These JET experimental observations are corroborated by modelling with the local (or 'flux-tube') version of the state-of-the-art gyrokinetic code Gene [START_REF] Jenko | Electron temperature gradient driven turbulence[END_REF] that successfully reproduced turbulence properties for various tokamak experiments (see e.g. Refs. [6,7,[START_REF] Citrin | Electromagnetic stabilization of tokamak microturbulence in a high-β regime[END_REF][START_REF] Görler | On the Validation of Gyrokinetic L-mode Simulations[END_REF][START_REF] Mazzi | Impact of fast ions on a trapped-electron-mode dominated plasma in a JT-60U hybrid scenario[END_REF]). For analysis of JET pulse #94701, we choose a flux-tube around ρ tor = 0.23 (R = 3.25 m) as this region is characterized by large temperature gradients and a significant population of fast ions and TAEs. Table 1 summarizes the modelling input data. An equivalent Maxwellian distribution was used to represent the population of fast deuterons. Linear stability analysis shows that, as expected, ITG modes dominate the spectrum, peaking at k y ρ s ≈ 0.45. Here, k y is the binormal ('poloidal') wavenumber and ρ s is the characteristic Larmor radius of thermal ions at the reference sound speed defined in this paper as c s = (T e /m p ) 1/2 with m p the proton mass.

Because of the uncertainties in the fast ion profiles calculated by TRANSP in the presence of destabilized Alfvén eigenmodes and limitations of the Maxwellian approximation to represent fully the detailed fast ion distribution, a scan over the normalized fast ion logarithmic pressure gradient, R/L p FD was performed. This allows to compensate the discrepancy induced by the derivative in velocity space by fine-tuning the radial gradient of the fast ion temperature and density profiles. Furthermore, as a part of this assessment, we find that for R/L p FD > 10.7 fast D ions destabilize ballooned modes with frequencies f ≈ 200 kHz and k y ρ s between 0.025 and 0.05, corresponding to n = 4 -7. As all these mode characteristics match well with those of the experimentally measured TAEs (Figure 2(e)), these high-frequency modes appearing in the Gene simulations are identified as fast-ion-driven TAEs.

This analysis was followed by nonlinear Gene simulations. Figure 3(a) shows the heat diffusivity attributed to electrostatic fluctuations for both thermal ion species, χ ES (D) and χ ES ( 3 He). Figures 3(b)-(d) illustrate the amplitudes of the Fourier components normalized to their maximum, φn of the perturbed electrostatic potential φ (averaged over the radial and parallel directions) as a function of frequency and k y ρ s . The frequency dependence of φn , additionally averaged over k y ρ s , is illustrated in Figures 3(e)-(g). The turbulence spectra reveal that for R/L p F D ≈ 6, the dominant modes are characterized by low frequencies and peak around k y ρ s ≈ 0.2. These k y values are somewhat lower than in linear simulations. This is typical for a local gyrokinetic modelling and linked to the inverse cascade of energy, in which modes with larger k y transfer energy to modes with lower k y . At larger fast ion gradients, e.g. at R/L p F D ≈ 16, the turbulence structure is substantially modified. The dominant modes are shifted towards higher frequencies f ≈ 200 kHz and lower k y ρ s ≈ 0.025 -0.05, matching closely the characteristics of the observed TAEs.

Figure 3(a) clearly shows a strong stabilizing effect of fast ions on ion-scale turbulence when fast ion driven TAEs appear in nonlinear Gene simulations; see, e.g., Figure 3(g) computed for R/L p F D ≈ 16. Under these conditions, the computed diffusivities of both thermal ion species are reduced by more than 95% compared to the values obtained without fast ions, and agree well with χ i computed by TRANSP (the black dashed line in Figure 3(a)). However, the predicted Gene value χ e ≈ 3.5 m 2 s -1 is significantly larger than the electron heat diffusivity from TRANSP (χ e ≈ 0.5 m 2 s -1 ). To obtain better agreement between Gene and power balance results, an additional set of non-linear simulations was performed, by varying the input parameters within experimental uncertainties. The consistency in both χ i and χ e was simultaneously achieved at R/L pF D ≈ 14 with temperature gradients reduced by 15%, see Figure 3(c) and (f), yielding χ e ≈ 0.7 m 2 s -1 . This result supports the observation of high T e0 in JET pulse #94701, in contrast to earlier NSTX results where the electron confinement was strongly degraded in the presence ofAlfvén eigenmodes [START_REF] Gorelenkov | Anomalous electron transport due to multiple high frequency beam ion driven Alfvén eigenmodes[END_REF].

We note that modes at f ≈ 20 -30 kHz were observed in #94701. Modes at similar frequencies also appeared in Gene simulations at R/L p F D ≈ 14 and 16. The bispectral analysis of the electrostatic potential indicates a nonlinear coupling between these modes and TAEs, similar to results reported in [START_REF] Crocker | Three-wave interactions between fast-ion modes in the national spherical torus experiment[END_REF]. However, since turbulence suppression is reached in Gene simulations already at R/L p F D ≈ 11, when such a mode coupling is not observed, we exclude these low-frequency-modes as a main cause for the observed confinement improvement.

Figures 4(a)-(c) illustrate the dependence of the perturbed electrostatic potential φ on the radial (x/ρ s ) and binormal (y/ρ s ) coordinates. In the absence of fast ions (Fig. 4(a)), small-scale turbulent eddies predominately elongated along the radial direction, typical for ITG-induced transport, are seen. In contrast, Figs. 4(b) and(c) show the sudden appearance of intense poloidally-oriented and radially-sheared zonal flows in GENE simulations at R/L p F D ≈ 14 and R/L p F D ≈ 16, corresponding to plasmas with MeV ions and destabilized TAEs. In these cases, the zonal flow shearing rate respectively reaches γ E×B, zonal ≈ 0.84 c s /a and 0.86 c s /a, almost twice the value without fast ions (0.48 c s /a). These zonal flows strongly affect the ITG nonlinear saturation [START_REF] Diamond | Zonal flows in plasmas -a review[END_REF] and produce a radial de-correlation of the turbulent eddies [START_REF] Biglari | Influence of sheared poloidal rotation on edge turbulence[END_REF], thereby leading to a reduction in the turbulent transport. Furthermore, Figs. 4 (e)-(f) illustrate the appearance of zonal structures also in the perturbed magnetic potential with unstable TAEs, consistent with earlier analytical studies [START_REF] Qiu | Effects of energetic particles on zonal flow generation by toroidal Alfvén eigenmode[END_REF]. The strong coupling between fast ion driven modes and zonal components, previously shown in simplified contexts [START_REF] Todo | Nonlinear magnetohydrodynamic effects on Alfvén eigenmode evolution and zonal flow generation[END_REF][START_REF] Qiu | Effects of energetic particles on zonal flow generation by toroidal Alfvén eigenmode[END_REF], is further highlighted by the bispectral analysis. This analysis suggests a net energy transfer from TAEs to the zonal modes, see Supplementary Figure S6.

The robustness of Gene results was further checked against the uncertainty in the input R/L Ti , arising from a lack of T i data at ρ tor < 0.2 (see also Methods). The mechanism is reproducible and robust against unavoidable assumptions and simplifications of the local approach.

Experimental evidence for minimized density fluctuations

Modeling predictions of the turbulence suppression are supported by the observation of reduced amplitude of the density fluctuations in these JET plasmas, in particular, in pulse #95669 performed at very similar plasma conditions as #94701. Figures S3(b) and (c) in the Supplementary Information show that much smaller density fluctuation amplitudes were observed in the plasma phase with MeV-range fast ions and fast ion driven TAEs, as compared to the NBI-only phase. As outlined earlier, partial turbulence reduction in the presence of moderately energetic ions (of a few hundred keV) does not necessarily require TAEs and has been extensively reported [6][7][8][START_REF] Di Siena | Electromagnetic turbulence suppression by energetic particle driven modes[END_REF]. This has also been observed in this series of JET experiments, e.g. in pulse #95672 (see Supplementary Fig. S3), where the application of the three-ion scheme with P ICRF = 2.1 MW resulted in monster sawteeth, while TAEs were not destabilized. However, a much stronger reduction was measured in the same JET pulse when ICRF power was increased to 6 MW such that MeV-range ions and fast ion driven TAEs were present in the plasma. Supplementary Figure S4(c) clearly illustrates the efficacy of the turbulence suppression mechanism discovered in these JET experiments.

A promising result for ITER and future fusion devices

These promising JET results indicate that unexpected favourable conditions might be realized in ITER and future fusion power plants, where alpha particles provide a strong source of core electron heating [START_REF] Stix | Heating of toroidal plasmas by neutral injection[END_REF]. As fast alphas can simultaneously destabilize TAEs at large plasma volumes in ITER [START_REF] Pinches | Energetic ions in ITER plasmas[END_REF], this complex mechanism -if present -would provide a significant increase in plasma confinement. However, TAEs can also enhance the fast particle transport, and thus the final result would depend on a delicate balance between different competing effects. Predicting the actual TAE mode stability and the efficiency of turbulence suppression in ITER plasmas is challenging because of the sensitivity to radial profiles of background quantities or the energetic particle distributions. Therefore, whether this favorable mechanism of turbulence stabilization will actually materialize in ITER plasmas is an open question. However, the results shown in this paper pave the way to an enhanced realization of fusion as an energy source, which is worth exploring further both experimentally and theoretically in the future.

Tables Table 1: Plasma parameters used in the Gene modelling of JET pulse #94701 at ρtor = 0.23 and t ≈ 9.6 s. Here, R0 is the major radius, B0 the on-axis magnetic field, the inverse aspect ratio, q and ŝ the local safety factor and magnetic shear, nj and Tj the local density and temperature of various plasma species (electrons, thermal D and 3 He ions, fast D ions), R/Ln,T the normalized logarithmic density and temperature gradient, βe the electron-beta, ν * the normalized collision frequency. The reported input parameters are common to all the numerical Gene simulation cases. The various cases, however, differ essentially in the fast ion pressure gradient R/Lp FD , whose values are displayed e.g. in Figure 3. In tokamaks, charged particles are confined by a combination of toroidal and poloidal magnetic fields that generate nested surfaces of constant magnetic flux. Passing particles (blue trajectory) circulate along magnetic field lines in the toroidal direction ϕ, while trapped particles (red trajectory) periodically reverse their parallel velocity with respect to the confining magnetic field B. For simplicity, the rapid gyromotion of charged particles is not shown and only the guiding center motion is displayed. In this paper, the local transport analysis focuses on the magnetic flux surface corresponding to ρtor = 0.23 (in green), where ρtor is the square root of the toroidal magnetic flux normalized to its value at the plasma boundary. The top left inset illustrates the field-aligned set of coordinates used in the flux-tube version of the Gene code, with (x, y, z) the radial, binormal and parallel coordinate, respectively. The top right inset shows a poloidal cross-section of the plasma in the (R,Z) plane for the reported JET experiments, together with the studied flux surface ρtor = 0.23 (green) and the plasma boundary ρtor = 1 (grey). Here, R is the radial distance from the torus centre and Z is the vertical coordinate. Note that ρtor = 0.23 corresponds to R = 3.25 m at the outboard side. (i) display the computed electron, χe and thermal ion, χi heat diffusivities. The error bands in these panels take into account uncertainties in the temperature profile data, different fittings to Ti and its extrapolation to the magnetic axis, including uncertainties in the power balance. The red dotted and dashed lines in (g) show the predicted Ti profiles for pulse #94701, using the two χi profiles from panel (i) as input.

R 0 (m) B 0 (T) q ŝ n e (m -

Figure Legends/Captions

Methods

GENE simulation model and input parameters

In its flux-tube version, the Gene code solves the nonlinear gyrokinetic Vlasov equations [1] coupled to Maxwell's equations on a field-aligned set of spatial coordinates (for a schematic representation, see the top left inset of Figure 1). This system of coordinates allows the exploitation of the strong anisotropy of the turbulent fluctuations in parallel and perpendicular directions to the background magnetic field. Whereas the parallel spatial direction z employs a finite-difference solution technique, the perpendicular direction is treated with spectral methods, and hence the radial x and binormal y coordinates in the Fourier domain are generally referred as k x and k y respectively. The two velocity dimensions employed in Gene are the parallel velocity v || and the magnetic moment µ. A detailed description and derivation of the model equations employed in the flux-tube version can be found in Refs. [2,3].

In the present work, both perpendicular and parallel magnetic field fluctuations are computed, and the collisions are retained. In Gene, Maxwellian distribution functions are employed for all the particle populations, including fast D ions. For the fast ion distribution, effective density and temperature were calculated from the TRANSP distribution. The number of grid points used in the nonlinear simulations is (n kx = 256, n ky = 48, n z = 32, n v || = 48, n µ = 64) for 4 different particle populations, i.e. electrons, thermal D and 3 He ions, and fast D ions. The minimum wavenumber in the binormal direction considered is k y,min ρ s = 0.025. The employed numerical discretization was chosen after extensive convergence tests, even in the nonlinear phase. It must be noted that in the simulations without fast particles, the velocity space has been relaxed to a less demanding numerical grid. Namely, the number of points in the µ-direction has been set to n µ = 16. Indeed, in order to accurately resolve the high-frequency resonances in fast ion velocity space, the non-equidistant Gauss-Legendre discretization in this direction had to be increased up to n µ = 64. Also, for a more comprehensive description of the numerical implementation techniques and adopted schemes employed in Gene, one can consult dedicated Ph.D. theses, e.g. [3], which are available on the Gene website. The input parameters that have been employed in the numerical analyses are reported in Table 1.

The energy flux observable is used to evaluate the radial energy transport computed by the Gene code and compare it with the experimental power balances. The definition of the flux-surface averaged energy flux in the radial direction is:

Q s = d 3 v 1 2 m s v 2 f s (x, v)v E×B (x)
where x = (x, y, z) and v = (v , µ), • denotes the flux-surface average,

f s (x, v) is the distribution function, v E×B = c B 2 0 B 0 × ∇ ξ is the generalized E × B drift velocity, with ξ = φ - v c
Ā + µ qs B the gyroaveraged modified potential (in which q s is the charge and the upper bar indicating a gyroaveraged quantity), and the subscript s refers to the species. From this definition, it is possible to appreciate the contributions to the energy flux of the electrostatic potential φ, the vector potential Ā and the parallel magnetic fluctuation B due to the v E×B drift. Therefore, it is possible to separate the electrostatic, which includes only the terms with φ, and the electromagnetic contributions, which includes the terms with Ā and B . In addition, it is also possible to determine the flux-surface averaged thermal diffusivities χ s from the energy flux as:

χ s = Q s n 0,s T 0,s ω T,s
where ω T,s = -a T0,s dTs dx is the normalized temperature gradient of the considered species s. For a more detailed derivation in Gene units of the energy flux and diffusivity, one can consult, e.g., [4].

Interpretive integrated modelling framework

Both pulses #94701 and #94704 were analysed through interpretive simulations performed with the TRANSP modelling suite [5] coupled with external heating modules NUBEAM (NBI) [6] and TORIC (ICRF) [7], and prepared with the OMFIT integrated modelling platform [8]. The interpretive analysis was based on the use of fitted profiles, including electron density and temperature. The fitting procedures for both quantities were based on high resolution Thomson scattering measurements, while the temperature was additionally constrained by electron cyclotron emission data, up to the availability of the experimental measurements. To extrapolate the ion temperature profiles to the magnetic axis, two different approaches were used for pulse #94701: i) assuming T i = T e in the region ρ tor ≤ 0.2 in combination with the available high-resolution T e data or ii) applying a global third-order polynomial fit in the range ρ tor 0.8 (with the additional constraint ∂T i (0)/∂r = 0). The latter fitting procedure is routinely employed in TRANSP modelling of JET experiments, with as a constraint the consistency between measured and computed quantities such as the neutron rate and the plasma stored energy. Both approaches for pulse #94701 provide very similar results. Indeed, for pulse #94704, only the third-order polynomial fitting procedure has been employed in the TRANSP simulation, with a good agreement (within ∼10%) with the measured neutron rate. The equilibrium used was an EFIT reconstruction constrained by magnetics and pressure profiles, i.e. including kinetic profiles as well as the contribution of fast ions.

Description of experimental measurements

Mirnov coils are used as a standard MHD diagnostic on almost all tokamak devices. The coils are installed within the vacuum vessel close to the plasma boundary and provide a measurement of the time derivative of the magnetic field. Magnetic spectrograms (Fourier decompositions of the Mirnov coil signal) can then be used to identify relevant oscillation frequencies associated with MHD activity. In JET a number of coil arrays with high frequency response are available, allowing activity in the Alfvén range to be observed. The radial localization of the modes was obtained using an X-mode reflectometer (see Supplementary information). The ion temperature profiles in this paper were obtained from Charge eXchange Recombination Spectroscopy (CXRS) measurements and electron temperature profiles from combined analysis of the electron cyclotron emission (ECE) and high resolution Thomson scattering (HRTS) diagnostics. The density profiles were taken from HRTS measurements, with the density normalized to match the line averaged density measured by a Far Infrared interferometer. The time resolved neutron yield in JET is measured using three fission chambers, containing 235 U and 238 U, located outside the vacuum vessel. The time-of-flight spectrometer for rate (TOFOR) measures the energy distribution of fusion-born neutrons by based on neutron time-of-flight measurements.

Figure 1 :

 1 Figure1: Schematic view of the tokamak geometry. In tokamaks, charged particles are confined by a combination of toroidal and poloidal magnetic fields that generate nested surfaces of constant magnetic flux. Passing particles (blue trajectory) circulate along magnetic field lines in the toroidal direction ϕ, while trapped particles (red trajectory) periodically reverse their parallel velocity with respect to the confining magnetic field B. For simplicity, the rapid gyromotion of charged particles is not shown and only the guiding center motion is displayed. In this paper, the local transport analysis focuses on the magnetic flux surface corresponding to ρtor = 0.23 (in green), where ρtor is the square root of the toroidal magnetic flux normalized to its value at the plasma boundary. The top left inset illustrates the field-aligned set of coordinates used in the flux-tube version of the Gene code, with (x, y, z) the radial, binormal and parallel coordinate, respectively. The top right inset shows a poloidal cross-section of the plasma in the (R,Z) plane for the reported JET experiments, together with the studied flux surface ρtor = 0.23 (green) and the plasma boundary ρtor = 1 (grey). Here, R is the radial distance from the torus centre and Z is the vertical coordinate. Note that ρtor = 0.23 corresponds to R = 3.25 m at the outboard side.

Figure 2 :

 2 Figure 2: Improved thermal ion confinement in fusion plasmas with MeV ions and unstable TAEs. This is illustrated by comparing two JET pulses #94704 (∼100 keV ions; blue lines) and #94701 (MeV-range ions; red lines) performed at the same operational conditions and the same total auxiliary heating power Paux = 14 MW, but differing in the characteristics of the fast ion population. (a) NBI and ICRF power. (b) Core electron temperature Te0. (c) Plasma stored energy Wp. (d) Neutron rate from D-D fusion reactions. (e) The magnetic coil spectrogram showing destabilized TAEs with different toroidal mode numbers, n in pulse #94701. Panels (f) and (g) show the measured radial profiles of the electron, Te and ion temperatures, Ti at Paux = 14 MW. The error bands in both panels represent a combination of the systematic diagnostic uncertainties and time-averaged (±0.1 s) statistical variations of the measurements. Panels (h) and(i) display the computed electron, χe and thermal ion, χi heat diffusivities. The error bands in these panels take into account uncertainties in the temperature profile data, different fittings to Ti and its extrapolation to the magnetic axis, including uncertainties in the power balance. The red dotted and dashed lines in (g) show the predicted Ti profiles for pulse #94701, using the two χi profiles from panel (i) as input.

Figure 3 :

 3 Figure 3: Thermal ion transport suppression in plasmas with MeV ions and fast ion driven modes. (a) The electrostatic thermal diffusivities of D and 3 He ions against the fast ion pressure gradient R/Lp FD as predicted by Gene.(b-d) Contour plots of the Fourier components normalized to their maximum, φn of the perturbed electrostatic potential φ (averaged over the radial and parallel directions) as a function of frequency and kyρs. (e-g) The frequency dependence of φn, additionally averaged over kyρs (indicated by the symbol ... ). In (f) and (g), the grey shaded areas indicate the experimentally measured TAE frequency range. In panel (a), the error bars of Gene time-dependent simulations represent the standard deviation obtained from the code calculation. The large error bars for the rightmost points can be attributed to the high-frequency and large-amplitude modulations of the energy fluxes due to strong fast ion driven fluctuations of the electrostatic potential. The horizontal dotted lines and the associated shaded areas in panel (a), represent thermal ion diffusivities and their standard deviations for the Gene modelling case without fast ions. The dashed line at the bottom of panel (a) shows the computed thermal ion heat diffusivity for JET pulse #94701, χi = 0.8 m 2 s -1 , as given by TRANSP, with the grey shaded area corresponding to the error band in Fig.2(i).

Figure 4 :

 4 Figure 4: Evidence for the nonlinear generation of zonal structures. Contour plots of the gyroaveraged perturbed electrostatic potential, φ, and the component of the vector potential parallel to the magnetic field, A || , for different fast ion pressure gradients R/Lp F D along the radial x and binormal y components. Panels (a) and (d) correspond to the modeling case without fast ions; panels (b) and (e) for the fast ion case with R/Lp F D ≈ 14; panels (c) and (f) for the fast ion case with R/Lp F D ≈ 16. All contour plots are normalized to the corresponding maximum value. This figure clearly shows the appearance of zonal structures leading to turbulence suppression in plasmas with MeV-range ions and fast ion driven TAEs.

Acknowledgements

We thank Matteo Baruzzo and Filomena Nave for the preparation and execution of JET experiments discussed in this paper; Elena de la Luna for her support in detailing the experimental diagnostics of JET; Aaron Ho for his assistance in processing the experimental data; Tobias Görler for providing an essential advice to ensure the correct numerical setup for the GENE modeling reported in this paper; Yann Camenen, Xavier Garbet and Andreas Bierwage for fruitful discussions about the gyrokinetic analyses; Gerardo Giruzzi for valuable suggestions on the article strategy. The simulations were performed on IRENE Joliot-Curie HPC system, in the framework of the PRACE projects IONFAST and AFIETC, and on CINECA Marconi HPC within the project GENE4EP. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under Grant agreement No 633053. The views and opinions express herein do not necessarily reflect those of the European Commission. Part of the work by Ye.O.K. and J.O. was also carried out in the framework of projects done for the ITER Scientist Fellow Network (ISFN).

Data availability

The JET experimental data is stored in the PPF (Processed Pulse File) system which is a centralised data storage and retrieval system for data derived from raw measurements within the JET Torus, and from other sources such as simulation programs. These data are fully available for the EUROfusion consortium members and can be accessed by non-members under request to EUROfusion.

Numerical data that support the outcome of this study are available from the corresponding author upon reasonable request.

Code availability

The research codes cited in the paper require a prior detailed knowledge of the implemented physics models and are under continuous development. The corresponding author can be contacted for any further information.

Author Contributions

The reported JET experiments were designed and coordinated by Ye.O.K., M.N., J.G. and J.O.; S.M. performed gyrokinetic simulations and subsequent analysis, with support of J.G., D.Z. and S.B.; input data for gyrokinetic modeling were provided by Ž. Š., G.S., Ye.O.K. and M.D.; Ž. Š. performed TRANSP modeling; J.G. performed power balance analysis and CRONOS simulations; M.D. provided analysis of the TAE radial location and the correlation reflectometer data; J.E. and A.S. provided neutron measurements data from TOFOR. The original manuscript was written by S.M., J.G., Ye.O.K., J.O. and D.Z. with feedback by all the authors. Major revisions of this manuscript were written by Ye.O.K., J.O., J.G. and S.M.

Competing Interests

The authors declare no competing interests.