Keywords: Semi-Lagrangian method, Runge-Kutta, plasma simulation, Vlasov equation méthode semi-Lagrangienne, Runge-Kutta, simulation de plasmas, équation de Vlasov

This work deals with the numerical solution of the Vlasov equation. This equation gives a kinetic description of the evolution of a plasma, and is coupled with Poisson's equation for the computation of the self-consistent electric field. The coupled model is non linear. A new semi-Lagrangian method, based on forward integration of the characteristics, is developed. The distribution function is updated on an eulerian grid, and the pseudo-particles located on the mesh's nodes follow the characteristics of the equation forward for one time step, and are deposited on the 16 nearest nodes. This is an explicit way of solving the Vlasov equation on a grid of the phase space, which makes it easier to develop high order time schemes than the backward method.

Introduction

Understanding the dynamics of charged particles in a plasma is of great importance for a large variety of physical phenomena, such as the confinement of strongly magnetized plasmas, or laser-plasma interaction problems for example. Thanks to recent developments in computational science and in numerical methods, meaningful comparisons between experience and numerics are becoming possible.

An accurate model for the motion of charged particles, is given by the Vlasov equation. It is based on a phase space description so that non-equilibrium dynamics can accurately be investigated. The unknown f (t, x, v) depends on the time t, the space x and the velocity v. The electromagnetic fields are computed self-consistently through the Maxwell or Poisson equations, which leads to the nonlinear Vlasov-Maxwell or Vlasov-Poisson system.

The numerical solution of such systems is most of the time performed using Particle In Cell (PIC) methods, in which the plasma is approximated by macro-particles (see [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF]). They are advanced in time with the electromagnetic fields which are computed on a grid. However, despite their capability to treat complex problems, PIC methods are inherently noisy, which becomes problematic when low density or highly turbulent regions are studied. Hence, numerical methods which discretize the Vlasov equation on a grid of the phase space can offer a good alternative to PIC methods (see [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF][START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF][START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF][START_REF] Carillo | Non oscillatory interpolation methods applied to Vlasov-based models[END_REF]). The so-called Eulerian methods can deal with strongly nonlinear processes without additional complexity, and are well suited for parallel computation (see [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF]). Moreover, semi-Lagrangian methods which have first been introduced in meteorology (see [START_REF] Staniforth | Semi-Lagrangian integration schemes for atmospheric models -A review[END_REF][START_REF] Zerroukat | A monotonic and positive-definite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF][START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF]), try to take advantage of both Lagrangian and Eulerian approaches. Indeed, they allow a relatively accurate description of the phase space using a fixed mesh and avoid traditional step size restriction using the invariance of the distribution function along the trajectories. Standard semi-Lagrangian methods calculate departure points of the characteristics ending at the grid point backward in time; an interpolation step enables to update the unknown.

In this work, we consider the numerical resolution of the two-dimensional Vlasov equation on a mesh of the phase space using a forward semi-Lagrangian numerical scheme. In the present method, the characteristics curves are advanced in time and a deposition procedure on the phase space grid, similar to the procedure used in PIC methods for the configuration space only, enables to update the distribution function.

One of the main cause for concern with semi-Lagrangian methods, is computational costs. With Backward semi-Lagrangian methods (BSL), the fields have to be computed iteratively, with Newton fixed point methods, or prediction correction algorithms. This is due to an implicit way of solving the characteristics (see [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF] for details). This strategy makes high order resolution quite difficult and expensive. Making the problem explicit enables to get rid of iterative methods for the characteristics, and to use for example high order Runge Kutta methods more easily. This is one of the main advantages of the present Forward semi-Lagrangian (FSL) method. Once the new position of the particles computed, a remapping (or a deposition) step has to be performed. This issue is achieved using cubic spline polynomials which deposit the contribution of the Lagrangian particles on the uniform Eulerian mesh. This step is similar to the deposition step which occurs in PIC codes but in our case, the deposition is performed in all the phase space grid. Similarities can also be found in strategies developed in [START_REF] Cotter | Reich The remapped particle-mesh semi-Lagrangian advection scheme[END_REF][START_REF] Nair | A forward-trajectory global semi-Lagrangian transport scheme[END_REF][START_REF] Reich | An explicit and conservative remapping strategy for semi-Lagrangian advection[END_REF] for meteorology applications.

In order to take benefit from the advantages of PIC and semi-Lagrangian methods, and since the two methods (PIC and FSL) really look like each other, except the deposition step, we have also developed a hybrid method, where the deposition step is not performed at each time step, but every T time steps. During the other time steps, the fields are computed directly at the new position of the particles. It shall be noticed that however the present method is not a real PIC method, since the particle weights are not constant. Indeed, in this method based on a description of the unknown using cubic spline polynomials, the spline coefficients play the role of the particle weights, and are updated at each time step. This kind of hybrid approach has been developed recently in a slightly different framework in [START_REF] Vadlamani | The particle-continuum method: an algorithmic unification particle in cell and continuum methods[END_REF] inspired by [START_REF] Denavit | Numerical simulation of plasmas with periodic smoothing in phase space[END_REF].

This paper is organized as follows. In the next section, the two Vlasov equations which will be dealt with are presented. Then, we shall introduce the Forward semi-Lagrangian (FSL) method, always regarding it comparatively to Backward semi-Lagrangian (BSL) methods. Afterward, numerical results for several test cases are shown and discussed. Eventually, some specific details are given in two appendices, one for the computation of an exact solution to the Landau damping problem, and the other about the solution of the Poisson equation for the Guiding-Center model.

Models in plasma physics

In this section, we briefly present two typical reduced models from plasma physics for the description of the time evolution of charged particles. These two-dimensional models are relevant for more complex problems we are interested in and shall be used to validate our new method.

Vlasov-Poisson model

We consider here the classical 1D Vlasov-Poisson model, the unknown of which f = f (t, x, v) is the electron distribution function. It depends on the space variable x ∈ [0, L] where L > 0 is the size of the domain, the velocity variable v ∈ IR and the time t ≥ 0. The Vlasov equation which translates the invariance of the distribution function along the characteristics then writes

∂f ∂t + v∂ x f + E(t, x)∂ v f = 0, (2.1)
with a given initial condition f (0,

x, v) = f 0 (x, v). The self-consistent electric field E(t, x) is computed thanks to the distribution function f ∂ x E(t, x) = I R f (t, x, v)dv -ρ i , L 0 E(t, x)dx = 0, (2.2)
where ρ i denotes the ion density which forms a uniform and motionless background in the plasma.

The Vlasov-Poisson model constitutes a nonlinear self-consistent system as the electric field determines f with (2.1) and is in turn determined by it in (2.2). It presents several conserved quantities as the total number of particles, the L p norms (p ≥ 1) defined by f L p = ( |f | p dxdv) 1/p , the momentum and the total energy, as follows:

d dt f (t, x, v)dxdv = d dt f (t) L p = d dt vf (t, x, v)dxdv = d dt v 2 f (t, x, v)dxdv + E(t, x) 2 dxdv = 0.
One of the main features of the present work is to develop accurate numerical methods which are able to preserve these exactly or approximately these conserved quantities for long times.

Guiding-center model

We are also interested in other kinds of Vlasov equations. For instance, in the guiding-center approximation. Charged particles in magnetized tokamak plasmas can be modeled by the density f = f (t, x, y) in the 2 dimensional poloidal plane by

∂f ∂t + E ⊥ (x, y) • ∇f = 0, (2.3)
coupled self-consistently to Poisson's equation for the electric field which derives from a potential Φ = Φ(x, y)

-∆Φ(t, x, y) = f (t, x, y), E(t, x, y) = -∇Φ(t, x, y).

(2.4)

In equation (2.3), the advection term E ⊥ = (E y , -E x ) depends on (x, y) and the time-splitting cannot be simply applied like in the Vlasov-Poisson case. Hence, this simple model appears to be interesting in order to test numerical methods.

The guiding-center model (2.3)-(2.4) also presents conserved quantities as the total number of particles and L 2 norm of f (energy) and E (enstrophy)

d dt f (t, x, y)dxdy = d f (t) 2 L 2 dt = d E(t) 2 L 2 dt = 0. (2.5)

Characteristic curves

We can re-write Vlasov equations in a more general context by introducing the characteristic curves dX dt = U (X(t), t).

(2.6)

Let us introduce X(t, x, s) as the solution of this dynamical system, at time t whose value is x at time s. These are called the characteristics of the equation. With X(t) a solution of (2.6), we obtain:

d dt (f (X(t), t)) = ∂f ∂t + dX dt • ∇ X f = ∂f ∂t + U (X(t), t) • ∇ X f = 0. (2.7)
which means that f is constant along the characteristics. Using these notations, it can be written

f (X(t; x, s), t) = f (X(s; x, s), s) = f (x, s)
for any times t and s, and any phase space coordinate x. This is the key property used to define semi-Lagrangian methods for the solution of a discrete problem.

3 The forward semi-Lagrangian method

In this section, we present the different stages of the forward semi-Lagrangian method (FSL) and try to emphasize the differences with the traditional backward semi-Lagrangian method (BSL).

General algorithm

Let us consider a grid of the studied space (possibly phase-space) with N x and N y the number of points in the x direction [0, L x ] and in the y direction [0, L y ]. We then define ∆x = L x /N x , ∆y = L y /N y , x i = i∆x, y j = j∆y, for i = 0, .., N x and j = 0, .., N y . One important point of the present method is the definition of the approximate distribution functions which are projected on a cubic B-splines basis:

f (t, x, y) = k,l ω n k,l S(x -X 1 (t; x k , y l , t n ))S(v -X 2 (t; x k , y l , t n )), (3.8)
where X(t; x k , y l , t n ) = (X 1 , X 2 )(t, x k , y l , t n ) corresponds to the solution of the characteristics at time t (of the two dimensional system (2.6)) whose value at time t n was the grid point (x k , y l ). The cubic B-spline S is defined as follows

6S(x) =    (2 -|x|) 3 if 1 ≤ |x| ≤ 2, 4 -6x 2 + 3|x| 3 if 0 ≤ |x| ≤ 1, 0 otherwise.
In the expression (3.8), the weight w n k,l is associated to the particle located at the grid point (x k , y l ) at time t n ; it corresponds to the coefficient of the cubic spline determined by the following interpolation conditions

f (t n , x i , y j ) = k,l ω n k,l S(x i -X 1 (t n ; x k , y l , t n ))S(y j -X 2 (t n ; x k , y l , t n )) = k,l ω n k,l S(x i -x k )S(y j -y l ).
Adding boundary conditions (for example the value of the normal derivative of f at the boundaries, we obtain a set of linear systems in each direction from which the weights ω n k,l can be computed as in [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF][START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. We can now express the full algorithm for the forward semi-Lagrangian method

Step 0: Initialize f 0 i,j = f 0 (x i , y j )

Step 1: Compute the cubic splines coefficients ω 0 k,l such that

f 0 i,j = k,l ω 0 k,l S(x i -x k )S(y j -y l ),
Step 2: Integrate (2.6) from t n to t n+1 , given as initial data the grid points X(t n ) = (x k , y l ) to get X(t; x k , y l , t n ) for t ∈ [t n , t n+1 ], assuming the advection velocity U is known. We shall explain in the sequel how it is computed for our typical examples.

Step 3: Project on the phase space grid using (3.8) with t = t n+1 to get f n+1 i,j = f n+1 (x i , y j )

Step 4: Compute the cubic spline coefficients ω n+1 k,l such that

f n+1 i,j = k,l ω n+1 k,l S(x i -x k )S(y j -y l ).
Go to Step 2 for the next time step.

FSL: An explicit solution of the characteristics

For BSL, especially for the solution of the characteristics, it is possible to choose algorithms based on two time steps with field estimations at intermediate times. Generally, you have to use a fixed-point algorithm, a Newton-Raphson method (see [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF]), a prediction correction one or also Taylor expansions (see [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF]) in order to find the foot of the characteristics. This step of the global algorithm costs a lot (see [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF]). It is no longer needed in FSL, where the starting point of the characteristics is known so that traditional methods to solve ODEs, like Runge-Kutta algorithms can be incorporated to achieve high order accuracy in time. Let us show the details of this explicit solution of the characteristics, in Vlasov-Poisson and Guiding-Center models.

In both cases, the dynamical system (2.6) has to be solved. With FSL, X(t n ), U (X(t n ), t n ) are known. You can choose your favorite way of solving this system on each time step, since the initial conditions are explicit. This leads to the knowledge of X(t n+1 ) and U (X(t n+1 ), t n+1 ) so that Step 2 of the previous global algorithm is completed.

As examples of forward solvers for the characteristics curves, the second-order Verlet algorithm, Runge-Kutta 2 and Runge-Kutta 4 will be proposed for Vlasov-Poisson, and, as Verlet cannot be applied, only Runge-Kutta 2 and 4 will be used for the Guiding-Center model.

For Vlasov-Poisson, we denote by

X(t n ) = (X 1 (t n ), X 2 (t n )) = (x n , v n
) the mesh of the phase space, and

U (X(t n ), t n ) = (v n , E(x n , t n ))

the advection velocity. The Verlet algorithm can be written

Step 1: ∀k, l, v

n+ 1 2 k,l -v n l = ∆t 2 E(x n k , t n ), Step 2: ∀k, l, x n+1 k,l -x n k = ∆t v n+1/2 k,l ,
Step 3: compute the electric field at time t n+1

deposition of the particles x n+1 k,l on the spatial grid x i for the density ρ: ρ

(x i , t n+1 ) = k,l ω n k,l S(x i - x n+1
k,l ), like in a PIC method. -solve the Poisson equation on the grid x i : E(x i , t n+1 ).

Step 4: ∀k, l, v n+1 k,l -v

n+ 1 2 k,l = ∆t 2 E(x n+1 k,l , t n+1
). A second or fourth order Runge-Kutta algorithms can also be used to solve the characteristics curves of the Vlasov-Poisson system forward in time. The fourth order Runge-Kutta algorithm needs to compute intermediate values in time of the density and the electric field. Let us detail the algorithm omitting the indices k, l for the sake of simplicity

Step 1:

k 1 = (v n , E(x n , t n )) = (k 1 (1), k 1 (2)),
Step 2: compute the electric field at intermediate time t 1 :

deposition of the particles on the spatial grid x i for the density ρ: ρ(

x i , t 1 ) = k,l ω n k,l S[x i - (x n k + ∆t/2 k 1 (1))].
solve the Poisson equation on the grid x i : E(x i , t 1 ).

Step 3: compute

k 2 = (v n + ∆t 2 k 1 (2), E(x n + ∆t 2 k 1 (1), t 1 )
Step 4: compute the electric field at intermediate time t 2 :

deposition of the particles on the spatial grid x i for the density ρ: ρ(

x i , t 2 ) = k,l ω n k,l S[x i - (x n k + ∆t/2 k 2 (1))].
solve the Poisson equation on the grid x i : E(x i , t 2 ).

Step 5: compute

k 3 = (v n + ∆t 2 k 2 (2), E(x n + ∆t 2 k 2 (1), t 2 )
Step 6: compute the electric field at intermediate time t 3 :

deposition of the particles on the spatial grid x i for the density ρ: ρ

(x i , t 3 ) = k,l ω n k,l S[x i - (x n k + ∆t k 3 (1))].
solve the Poisson equation on the grid x i : E(x i , t 3 ).

Step 7: compute

k 4 = (v n + ∆t k 3 (2), E(x n + ∆t k 3 (1), t 3 ) Step 8: X n+1 -X n = ∆t 6 [k 1 + 2k 2 + 2k 3 + k 4 ]
In both Verlet and Runge-Kutta algorithms, the value of E at intermediate time steps is needed (step 3 for Verlet and steps 3, 5 and 7 for Runge-Kutta 4). This is achieved as in PIC algorithms by advancing the particles (which coincide at time t n with the mesh in this method) up to the required intermediate time.

Using a deposition step, the density is computed thanks to cubic splines of coefficients w n i on the mesh at the right time, and thus the electric field can also be computed at the same time thanks to the Poisson equation. Using an interpolation operator, the electric field is then evaluated at the required location (in steps 3, 5 and 7). Let us remark that this step involves a high order interpolation operator (cubic spline for example) which has been proved in our experiments to be more accurate than a linear interpolation (see section 4).

For the Guiding-Center equation, the explicit Euler method, and also Runge-Kutta type methods (of order 2, 3 and 4) have been implemented. There is no technical difficulty with computing high order methods. This is one of the general interests of forward methods. The time algorithm for solving the characteristics at the fourth order is similar to those presented in the Vlasov-Poisson case. However, there is a additional difficulty in the deposition step which enables to evaluate the density at intermediate time steps; indeed, the deposition is two-dimensional since the unknown does not depend on the velocity variable in this case.

Let us summarize the main steps of the second order Runge-Kutta method applied to the guiding center model of variables X n = (x n , y n ) and of advection field

U (X n , t n ) = E ⊥ (X n , t n ) Step 1: Xn+1 -X n = ∆tE ⊥ (X n , t n )
Step 2: Compute the electric field at time t n+1 -two-dimensional deposition of the particles on the spatial grid (x j , y i ) for the density ρ: ρ(x j , y i , t n+1 ) =

k ω n k S[x j -x n+1 k,l ]S[y i -y n+1 k,l
] -solve the two-dimensional Poisson equation on the grid x j : E(x j , y i , t n+1 ).

Step 3:

X n+1 -X n = ∆t 2 E ⊥ (X n , t n ) + E ⊥ ( Xn+1 , t n+1 )
Here, the numerical solution of the two-dimensional Poisson's equation is based on Fourier transform coupled with finite difference method. See details in Appendix II.

FSL -BSL Cubic Spline Interpolation

We are going to compare how spline coefficients are computed recurrently, for one dimensional transport problems, for the sake of simplicity.

FSL: deposition principle On our mesh, the grid points x i = i∆x, i = 0, .., N x at a time n can be regarded as particles. We have a distribution function which is projected onto a cubic splines basis. Thus, we know f (t n , x), ∀x, then the particles move forward, and we have to compute f (t n+1 , x i ), i = 0, ..., N x , reminding that f is constant along the characteristics, and that the particles follow characteristics between t n and t n+1 .

In fact, to each mesh point x i , a spline coefficient ω k is linked. The thing to understand, is that these coefficients are transported up to the deposition phase. The key is then to compute them recurrently as follows:

Deposition step

f n+1 (x i ) = k ω n k S(x i -X(t n+1 ; x k , t n )) = k/X(t n+1 ;x k ,t n )∈[xi-1,xi+2] ω n k S(x i -X(t n+1 ; x k , t n )),
Update of the splines coefficients ω n+1 k using the interpolation conditions

f n+1 (x i ) = i+2 k=i-1 ω n+1 k S(x i -x k ),
The number of points which actually take part in the new value of f n+1 (x i ) (here 4) is directly linked with the spline degree you choose. A p-Spline for example has a (p + 1) points support.

In a 1D way of regarding the problem, you can easily prove the mass conservation:

m n+1 = ∆x i f n+1 (x i ) = ∆x i k ω n k S(x i -X(t n+1 ; x k , t n )) = ∆x k ω n k = ∆x i f n (x i ) = m n .
Merely with the spline property of unit partition i S(x -x i ) = 1 for all x.

BSL: interpolation principle Let us introduce some notations. The foot of the characteristics X(t n , x i , t n+1 ) belongs to the interval [x l , x l+1 [. Then the reconstructed distribution function can be written Interpolation step

f n+1 (x i ) = f n (X(t n ; x i , t n+1 )) = l+2 k=l-1 ω n k S(X(t n ; x i , t n+1 ) -x k )
Update of the splines coefficients ω n+1 k using the interpolation conditions

f n+1 (x i ) = i+2 k=i-1 ω n+1 k S(x i -x k )
Figure 1: Principle of FSL (left) and BSL (right) for linear splines.

Here, we denoted by X(t n ; x i , t n+1 ) the foot of the characteristic coming from x i . The reader is referred to [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF][START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF] for more details on BSL interpolation.

In both cases, a linear system has to be solved, of equivalent complexity, so our method is as efficient at this level as BSL is.

Basic differences between FSL and BSL

Let us now explain the basic differences between forward and backward semi-Lagrangian methods. In both cases, a finite set of mesh points (x m ) m=1..N is used, and the values of the function f at the mesh points at a given time step t n are considered. The aim is to find the new values of f on the grid at the next time step t n+1 . BSL For BSL, in order to find the (n+ 1)-th value of f at x m , we follow the characteristic curve which goes through x m , backward in time, until time t n . The arrival point will be called the foot of the characteristics and does not necessarily coincide with a mesh point. Hence, we use any interpolation technique to compute f at this point, knowing all the values of the mesh at this time. This leads to the new value of f (x m ). Let us summarize: find the foot of the characteristics X(t n ) knowing X(t n+1 ) = x m (mesh point) interpolate using the grid function which is known at time t n . FSL For FSL, the principle is quite different. The characteristics beginning at time t n on the grid points are followed, during one time step, and the end of the characteristics (i.e. at time t n+1 ) is found. At this moment, the known value is deposited to the nearest grid points (depending on the chosen method). This deposition step is also performed in PIC codes on the spatial grid only, in order to get the sources for the computation of the electromagnetic field. Once every grid points has been followed, the new value of f is obtained by summing all contributions. The FSL method can be summarized as follows find the end of the characteristics X(t n+1 ) leaving from X(t n ) = x m (mesh point) deposit on the grid and compute the new particle weights.

Numerical results

This section is devoted to the numerical implementation of the forward semi-Lagrangian method. In particular, comparisons with the backward semi-Lagrangian method will be performed to validate the new approach.

Hill's equation

In order to check that high orders are really reached, a particularly easy model can be used, in which there are no self-consistent fields. This leads to a 1D model with an external force field written -a(t)x, where a is a given periodical function. The Vlasov equation becomes:

∂f ∂t + v∂ x f -a(t)x∂ v f = 0, (4.9)
The solution of this equation is seen through its characteristics, solutions of

dX dt = V, dV dt = -a(t)X (4.10)
thus, X is solution of Hill's equation:

d 2 X dt 2 + a(t)X = 0 (4.11)
Let's note that this equation can be written in a general way du dt = A(t)u, where A is a matrix valued periodic function. Since this is a 2D linear system, its solution is a 2D vector space and it is sufficient to find two independent solutions.

Let ω, ψ ∈ C 2 (R + , R), with ω(t) > 0 ∀t ∈ R + , so that ω is solution of the differential equation

d 2 ω dt 2 + a(t)ω - 1 ω 3 = 0 dψ dt = 1 ω 2 (4.12)
So u(t) = ω(t)e iψ(t) and v(t) = ω(t)e -iψ(t) are two independent solutions of Hill's equation (see [START_REF] Magnus | Hill's equation[END_REF] for more details) which can be determined numerically.

For this test case, the initial distribution function will be:

f 0 (x, v) = e -x 2 2ω 2 -ω 2 v 2 2 , ∀(x, v) ∈ [-12, 12] 2 .
The associated solution f (x, v, t) will depend only on Aω(t). In particular, f will have the same periodicity as a and ω. This is what will be used for testing the code. For different orders (2 and 4), and different ∆t, x rms (t) = x 2 f (x, v, t)dxdv will be displayed on Fig 2. This function should be periodic, and thus should reach the same test value x rms (0) at each period. The error will be measured between the ten first computed values and the exact one, for the ten first periods. Then these errors will be summed, so that a L 1 norm of the error is dealt with:

err = k=10 k=0 e k , with e k = |x rms (2kπ) -x rms (0)|.
The order of the method is checked in Figure 2. Note that N x = N v = 1024 to make sure that convergence is achieved for the interpolation step. The expected order is achieved for a certain ∆t interval. If ∆t becomes too small, a kind of saturation happens. This is due to the term in h m+1 ∆t (where h = ∆x = ∆v) in the theoretical estimation of the error for backward methods ( [START_REF] Besse | Convergence of classes of high order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF]), which becomes too high and prevents us from keeping the correct order. A forthcoming paper will try to do the same kind of error estimation for the forward method.

Vlasov-Poisson case

Landau damping The initial condition associated to the scaled Vlasov-Poisson equation (2.1)-(2.2) has the following form where k = 0.5 is the wave number and α = 0.001 is the amplitude of the perturbation, so that linear regimes are considered here. A cartesian mesh is used to represent the phase space with a computational domain [0, 2π/k] × [-v max , v max ], v max = 6. The number of mesh points in the spatial and velocity directions is designated by N x = 64 and N v = 64 respectively. Finally, the time step is equal to ∆t = 0.1 and the Verlet algorithm is used to compute the characteristics.

f 0 (x, v) = 1 √ 2π exp(-v 2 /2)(1 + α cos(kx)), (x, v) ∈ [0, 2π/k] × IR, ( 4
In this context, it is possible to find an exact value of the dominant mode solution of the linearized Vlasov-Poisson equation (see Appendix I for some details). The exact electric field corresponding to the dominant mode reads E(x, t) = 4α × 0.3677e -0.1533t sin(0.5x) cos(1.4156t -0.5326245).

On Fig. 3, the analytical solution of the L 2 norm of the electric field and the implemented one are plotted. It can be observed that the two curves are very close to each other. In particular, the damping rate and the frequency of the wave are well recovered (γ = -0.1533 and ω = 1.4156) by the method. Similar precision is achieved for different values of k leading to different value of the damping rate and of the frequency (see Fig. 3).

The recurrence effect that occurs with the present velocity discretization on a uniform grid, at T R ≈ 80 ω -1 p can also be remarked. This value is in good agreement with the theoretical recurrence time which can be predicted in the free-streaming case (see [START_REF] Nakamura | Cubic interpolated propagation scheme for solving the hyper-dimensional Vlasov-Poisson equation in phase space[END_REF]) T R = 2π k∆v . This test case has also been solved with the "hybrid" method in which the deposition step is only performed every T time steps. In all other steps, the remapping (or deposition) step is not performed, therefore, it can be linked with a PIC method. As it was already said, it is not really a PIC method since the spline coefficients are different on the phase space grid and are updated at each remapping step, whereas in classical PIC methods, these coefficients (called weights) are constant equal to n0

Npart where N part is the number of particles. On Fig. 4, the electric field is plotted again, for different values of T , and ∆t = 0.1, with N x = N v = 128 points. As expected, the method works well, even for large values of T . Only a kind of saturation can be observed, and it can be seen that values smaller than 2 -18 are not well treated, because of the lack of accuracy of the hybrid method. Nevertheless, the results are convincing: the computation gets faster as T gets larger, and a good accuracy is still reached. Two stream instability This test case simulates two beams with opposit velocities that encounter (see [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Nakamura | Cubic interpolated propagation scheme for solving the hyper-dimensional Vlasov-Poisson equation in phase space[END_REF]). The corresponding initial condition can be given by

f 0 (x, v) = M(v)v 2 [1 -α cos(kx)],
with k = 0.5 and α = 0.05. The computational domain is [0, 2π/k] × [-9, 9] which is sampled by N x = N v = 128 points. The Verlet algorithm is used to compute the characteristics with ∆t = 0.5. We are interested in the following diagnostics: the first three modes of the electric field, the electric energy 1/2 E(t) 2 L 2 and the time evolution of the phase space distribution function. On Fig. 5, we plot the time history of the first three Fourier modes of the electric field:

|E 1 |, |E 2 |, |E 3 |
denotes the amplitudes of Ê(k = 0.5) , Ê(k = 1) and Ê(k = 1.5) respectively. We observe that after an initial phase, the first mode exponentially increases to reach its maximum at T ≈ 18 ω -1 p . After this phase and until the end of the simulation, a periodic behavior is observed which translates the oscillation of the trapped particles in the electric field; in particular, a vortex rotates with a period of about 18 ω -1 p . The other modes |E 2 | and |E 3 | also grow exponentially and oscillate after the saturation. However, their amplitude remains inferior to that of the first mode. Similar observations can be performed for the electric energy which reaches its maximum at T ≈ 18 ω -1 p after an important and fast increase (from t = 8 to t = 18 ω -1 p ). This test case was also solved with the hybrid method to test the capability in the nonlinear regime. On Fig. 6, the first Fourier mode of the electric field is displayed for different T , with ∆t = 0.5, and 128 points in each direction. It can be observed that during the first phase, which is a linear one, even for big T , the results are quite accurate for all values of T . The hybrid method seems to have more difficulty after this linear phase. Numerical noise, one of the main drawbacks of PIC methods can be observed as T gets higher. The phenomenon can be understood looking at the distribution function. On Fig. 6 the noise clearly appears on f . Noisy values quickly reach high values which prevent the method from being accurate enough. They are more important in our hybrid method than in classical PIC ones, because of the deposition step, where the particles weight play their role. Indeed, if particles with very different weights are located at the same place, the deposition does not take into account properly the particles of low weights compared to those of heavy weights. On Fig. 6, it can be seen that the vortex which appears at the middle of the distribution and should stay there slowly leaves out of the domain. This can be explained as a kind of diffusion. Nevertheless, if T remains very little (2 -4), the results are really good. It can also be noticed that ∆t plays an important role, actually, when a smaller ∆t is chosen, the results remain good for bigger T . As an example, if ∆t = 0.1, results remain acceptable until T = 16. - Bump on tail Next, we can apply the scheme to the bump-on-tail instability test case for which the initial condition writes

f 0 (x, v) = f (v)[1 + α cos(kx)], with f (v) = n p exp(-v 2 /2) + n b exp - |v -u| 2 2v 2 t
on the interval [0, 20π], with periodic conditions in space. The initial condition f 0 is a Maxwellian distribution function which has a bump on the Maxwell distribution tail; the parameters of this bump are the following

n p = 9 10(2π) 1/2 , n b = 2 10(2π) 1/2 , u = 4.5, v t = 0.5,
whereas the numerical parameters are N x = 128, N v = 128, v max = 9, ∆t = 0.5. The Runge-Kutta 4 algorithm is used to compute the characteristics.

We are interested in the time evolution of the spatially integrated distribution function

F (t, v) = 20π 0 f (t, x, v)dx,
and in the time history of the electric energy 1/2 E(t) 2 L 2 . For this latter diagnostic, we expect oscillatory behavior of period equal to 1.05; moreover, since an instability will be declared, the electric energy has to increase up to saturation at t ≈ 20.95 and to converge for large times to 36% of its highest value (see [START_REF] Nakamura | Cubic interpolated propagation scheme for solving the hyper-dimensional Vlasov-Poisson equation in phase space[END_REF][START_REF] Shoucri | A two-level implicit scheme for the numerical solution of the linearized vorticity equation[END_REF]).

On Figures 7 and8, we plot the electric energy as a function of time. We can observe that oscillations appear, the period of which can be evaluated to 1.; then the maximum value is reached at t ≈ 21 and the corresponding amplitude is about 9, which is in very good agreement with the results presented in [START_REF] Nakamura | Cubic interpolated propagation scheme for solving the hyper-dimensional Vlasov-Poisson equation in phase space[END_REF]. Then the amplitude of the electric energy decreases and presents a slower oscillation due to the particle trapping. Finally, it converges to an amplitude of about 2.8 which is very close to the predicted value. However, for very large times (at t = 250 ω -1 p ), FSL using lower time integrator algorithms (Runge-Kutta 2 or Verlet algorithms) leads to bad results (see Fig. 8). A very precise computation of the characteristics is required for this test and the use of Runge-Kutta 4 is crucial. The results obtained by BSL (see Fig. 7) are very close to those obtained by FSL using Runge-Kutta 4. Moreover, if ∆t increases (up to ∆t ≈ 0.75), BSL presents some unstable results whereas FSL remains stable up to ∆t ≈ 1. The use of high order time integrators - leads to an increase of accuracy but also to more stable results as ∆t increases. It has been remarked that linear interpolation of the electric field is not sufficient to obtain accurate results with FSL-Runge-Kutta 4 and cubic spline are used to that purpose. Fig. 9 and Fig. 10 shows the time development of the spatially integrated distribution function for FSL and BSL. We observe that very fast, the bump begins to be merged by the Maxwellian and a plateau is then formed at t ≈ 30 -40 ω -1 p .

Guiding-center case

Kelvin-Helmoltz instability In order to validate our guiding-center code, we used two test cases introduced in [START_REF] Shoucri | A two-level implicit scheme for the numerical solution of the linearized vorticity equation[END_REF] and [START_REF] Ghizzo | A Hilbert-Vlasov code for the study of high-frequency plasma beatwave accelerator[END_REF]. The instability is created choosing an appropriate ρ 1 which will perturb the solution around the equilibrium one (ρ 0 , φ 0 ). Using the the work of Shoucri, we will take:

ρ(x, y, t = 0) = sin(y) + 0.015 sin( y 2 ) cos(kx)

where k = 2π Lx and L x the length of the domain in the x-direction The numerical parameters are:

N x = N y = 128, ∆t = 0.5.
The domain size has an impact on the solution. The interval [0, 2π] will be used on the y-direction, and respectively L x = 7 and L x = 10 . This leads to real different configurations: With L x = 7, Shoucri proved that the stable case should be dealt with. That is what was observed with this code.

With L x = 10, the unstable case is faced. The results prove it on figure Fig. 11 and12. For this test case, the evolution of the energy E 2 dxdy and enstrophy ρ 2 dxdy will also be plotted on Fig. 13. These should be theoretically invariants of the system. Like for other semi-Lagrangian methods, the energy lowers during the first phase, which is the smoothing one, where micro-structures can not be solved properly. Nevertheless, the energy is well conserved. Moreover, on Fig. 13, FSL using second and fourth order Runge-Kutta's methods are compared to the BSL method. As observed in the bump-on-tail test, the Runge-Kutta 4 leads to more accurate results and is then very close to BSL. However, BSL seems to present slightly better behavior compared to FSL-Runge-Kutta 4. But FSL-Runge-Kutta 4 enables to simulate such complex problems using higher values of ∆t (see Fig. 14). We observed for example that the use of ∆t = 1 gives rise to very reasonable results since the L 2 norm of the electric field E decreases of about 4 %. Let us remark that BSL becomes unstable for ∆t ≥ 0.7.

Conclusion and perspectives

In this paper, we introduced the forward semi-Lagrangian method for Vlasov equations. The method has been tested on two different models, the one-dimensional Vlasov-Poisson one, and the guiding-center one. Different test cases have been simulated, and they are quite satisfying. The results are in some cases a bit less accurate, with respect to the conservation of invariants, than with the classical BSL method, but enables the use of very large time steps without being unstable and recovering all the expected aims. No iterative methods anymore, and high order time schemes can be use in a straightforward manner. The next step will be to test the method with the Vlasov-Maxwell model, in which we will try to solve properly the charge conservation problem, which is the ultimate goal. We will try to use PIC results about that conservation, for example in [START_REF] Barthelmé | Le problème de conservation de la charge dans le couplage des équations de Vlasov et de Maxwell[END_REF]. We will also try to prove theoretically the convergence of this method.

Appendix I: Linearized Vlasov Poisson and Landau damping

In classical plasma physics textbooks, only the dispersion relations are computed for the linearized Vlasov-Poisson equation. However using the Fourier and Laplace transforms as for the computation of the dispersion relation and inverting them, it is straightforward to obtain an exact expression for each mode of the electric field (and also the distribution function if needed). Note that each mode corresponds to a zero of the dispersion relation.

The solution of the Landau damping problem is obtained by solving the linearized Vlasov-Poisson equation with a perturbation around a Maxwellian equilibrium, which corresponds to the initial condition f 0 (x, v) = (1 + ǫ cos(kx))/ √ 2πe -v 2 /2 . Let us introduce the plasma dispersion function Z of Fried and Conte [START_REF] Fried | The plasma dispersion function[END_REF]

Z(η) = √ πe -η 2 [i -erf i(η)], where erf i(η) = 2 π η 0 e t 2 dt.
We also have, Z ′ (η) = -2(ηZ(η) + 1). Then, denoting by E(k, t) the Fourier transform of E and by E(k, ω) the Laplace transform of E, the electric field, solution of linearized Vlasov-Poisson satisfies:

E(k, t) = j Res ω=ωj E(k, ω)e -iωt
where

E(k, ω) = N (k, ω) D(k, ω) D(k, ω) = 1 - 1 2k 2 Z ′ ( ω √ 2k ), N (k, ω) = i 2 √ 2k 2 Z( ω √ 2k )
The dispersion relation corresponds to D(k, ω) = 0.

For each fixed k, this equation has different roots ω j , and to which are associated the residues defining E(k, t) that can be computed with Maple. These residues that in fact the values

N (k, ω j ) ∂D ∂ω (k, ω) . ( 6.14)
Let us denote by ω r = Re(ω j ), ω i = Im(ω j ), r will be the amplitude of (6.14) and ϕ its phase.

Remark: For each root ω = ω r + iω i , linked to re iϕ , there is another: -ω r + iω i linked to re -iϕ . Then keeping only the roots in which ω i is the largest, which are the dominating ones after a short time, we get: E(k, t) ≈ re iϕ e -i(ωr +iωi)t + re -iϕ e -i(-ωr +iωi)t = 2re ωit cos(ω r t -ϕ)

Taking the inverse Fourier transform, we finally get an analytical expression for the dominating mode of the electric field, which we use to benchmark our numerical solution: E(x, t) ≈ 4ǫre ωit sin(kx) cos(ω r t -ϕ) Remark This is not the exact solution, because we have kept only the highest Laplace mode. Nevertheless, after about one period in time, this is an excellent approximation of E, because the other modes decay very fast. We have to solve:

-∆φ(x, y) = ρ(x, y)

We use a Fourier transform in the x direction. This leads, for i ∈ [1, N x ]:

∂ 2 φ i (y) ∂y 2 = ξ 2 φ i (y) + ρ i (y)
Let us introduce a notation. Thanks to Taylor Young formula, we have: This is nothing but the solution of a linear system

δu i = u i+1 -2u i + u i
A φ = R,
where A is a tridiagonal and symmetric matrix and R is a modified right hand side which allows to achieve a fourth order approximation.

Find E

To compute the electric field from the electric potential, we have to solve E = -∇φ. To achieve this task, a quadrature formula is used. In the x direction, which is the periodic one, a third order Simpson method is used where the (φ i ) i is given by previous step. There is no problem with extreme values, since the system is periodic. We then find the values of the electric field E solving another tridiagonal linear system. Whereas on the y-direction, Dirichlet conditions are imposed at the boundary. So, we can use the same strategy within the domain, but not on the two boundary points. Since we have a third order solution everywhere, we want to have the same order there, therefore we cannot be satisfied with a midpoint quadrature rule which is of second order. So we will add corrective terms, in order to gain one order accuracy. Here is how we do this. We want to find the precision of this method, thus, we would like to evaluate the following difference which is denoted by A: A = φ(x, 0) -φ(x, 1) -dy 2 (E(x, 0) + E(x, 1)) + dy

Figure 2 :

 2 Figure 2: Error as a function of ∆t for RK2 and RK4 (left) and x rms as a function of time, for ∆t = 2π/25, RK2 and RK4 (right).

Figure 3 :

 3 Figure 3: Linear Landau damping for k = 0.5 (left) and k = 0.4 (right)

Figure 4 :

 4 Figure 4: Linear Landau damping for k = 0.5 and for different number of T : from up to down and left to right: T = 1, T = 2, T = 16, T = 256.

Figure 5 :

 5 Figure 5: Two stream instability: Time evolution of the three first modes of the electric field (left) and of the electric energy (right).

Figure 6 :Figure 7 :

 67 Figure 6: Two stream instability: Time evolution of the first mode of the electric field (up and left) and distribution function at time t = 100 ω -1 p for T = 1, 4, 8.

- 1 pFigure 8 :

 18 Figure 8: Bump on tail instability: time evolution of the electric energy for FSL RK2 (left) and FSL Verlet (right).

Figure 9 :

 9 Figure 9: Bump on tail instability: time development of the spatially integrated distribution function for FSL (left) and BSL (right).

Figure 10 :

 10 Figure 10: Bump on tail instability: time development of the spatially integrated distribution function for FSL (left) and BSL (right).

Figure 11 : 1 p 4 . 3 . 1

 111431 Figure 11: Kelvin Helmholtz instability 1: distribution function at time t = 0, 30 ω -1 p

Figure 12 : 2 Figure 13 :

 12213 Figure 12: Kelvin Helmholtz instability 1: distribution function at time t = 50, 500 ω -1 p

Figure 14 :

 14 Figure 14: Kelvin Helmholtz instability 1: time history of L 2 norms of E (left) and of ρ (right). Comparison of the results for RK4 for different values of ∆t.

7

  Appendix II: Solution of Poisson in the Guided Center model 7.1 Find φ

xi+1 xi- 1 E

 1 (x, y)dx = -φ(x i+1 , y) + φ(x i-1 , y)

1 y ∆ 2 y 12 ξ 2 ( 2 y 12 [

 12122212 , y)dy = -φ(x, 1) + φ(x, 0) ≈ dy 2 (E(x, 0) + E(x, 1)) -dy 2 12 (ρ(x, 1) -ρ(x, 0)) + f φ where f φ = F -φ(ξ, 1) -φ(ξ, 0)) = ∆ -∂ xx (φ(x, 1) -φ(x, 0))] .

1 ) 6 ∂ 2 ∂y 2 E∂y 2 E

 16222 -ρ(x, 0)) -f φ .Using the Poisson equation and Taylor expansion, we haveA + E(x, 0) + E(x, 1) = -2 dy (φ(x, 1) -φ(x, 0)) + dy 6 (ρ(x, 1) -ρ(x, 0)) + dy 6 (∂ xx (φ(x, 1) -φ(x, 0)) = -2 dy (φ(x, 1) -φ(x, 0)) -dy 6 (∂ yy (φ(x, 1) -φ(x, 0)) = -2 dy (φ(x, 1) -φ(x, 0)) -dy 2 6 (∂ yyy φ(x, ξ 1 ) + O(∆ 3 y )where ξ 1 ∈ [y 0 , y 1 ]. Thus, we finally obtainA + E(x, 0) + E(x, 1) = -2 dy (φ(x, 1) -φ(x, 0)) + dy 2 (x, ξ 1 ) + O(∆ 3 y )Moreover, classical quadrature theory gives us the existence of ξ 2 ∈ [y 0 , y 1 ] such as (x, ξ 2 ).Replacing E(x, ξ 1 ) by E(x, ξ 2 ), which is of first order in our computation leads toA + E(x, 0) + E(x, 1) = -2 dy (φ(x, 1) -φ(x, 0)) -2dy y1 y0 E(x, y)dy + (E(x, 0) + E(x, 1)) + O(∆ 3 y ) so that A = O(∆ 3 y ) which is what was expected.

  (ξ 2 φ i + ρ i ) + O(∆ 2 y ) = ξ 2 φ i + ρ i + ∆y 2 12 (ξ 2 δ φ i + δ ρ i ) + O(∆ 2 y ).

				∆y 2	-1	= 1 +	∆y 2 12	∂ 2 ∂y 2	∂ 2 u i ∂y 2 + O(∆ 2 y ).
	Let us apply this to φ i
	∆y 2 12 ∂y 2 Now, factorizing all of it, we get: ∂ 2 δ φ i = 1 +
	φ i+1 1 -	ξ 2 ∆ 2 y 12	+ φ i -2 +	10ξ 2 ∆ 2 y 12	+ φ i-1 1 -	ξ 2 ∆ 2

y 12 = ∆y 2 ( ρ i+1 + 10 ρ i + ρ i-1 ) + O(∆ 4 y ).

RR n 6727


S Mazzi

J Garcia

D Zarzoso

Ye O Kazakov

J Ongena

M Nocente

M Dreval

Ž Štancar

G Szepesi

J Eriksson

A Sahlberg

S Benkadda

Towards enhanced performance in fusion plasmas via turbulence suppression by MeV ions

Megaelectron volt (MeV) alpha particles will be the main source of plasma heating in magnetic confinement fusion reactors. Yet, instead of heating fuel ions, most of the energy of alpha particles is transferred to electrons. Furthermore, alpha particles can also excite Alfvénic instabilities, previously considered as detrimental. Contrary to expectations, we demonstrate efficient ion heating in the presence of MeV ions and strong fast-ion driven Alfvénic instabilities in recent experiments on the Joint European Torus (JET). Detailed transport analysis of these experiments with state-of-the-art modeling tools explains the observations. Here we show a novel type of turbulence suppression and improved energy insulation in plasmas with MeV ions and fully developed Alfvénic activities through a complex multi-scale mechanism that generates large-scale zonal flows. This mechanism 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 [START_REF] Ongena | Magnetic-confinement fusion[END_REF]. The largest fusion device currently in operation is the Joint European Torus (JET) [START_REF] Litaudon | Overview of the JET results in support to ITER[END_REF]. An even larger experimental device ITER has started its assembly phase in July 2020, aiming to demonstrate the scientific and technological matureness of magnetic confinement fusion [START_REF] Shimada | Chapter 1: Overview and summary[END_REF].

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 to avoid damage of the wall of the device. These gradients create instabilities and turbulence, often with very different spatio-temporal scale lengths [START_REF] Doyle | Chapter 2: Plasma confinement and transport[END_REF], that result in energy and particle transport. 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 [START_REF] Romanelli | Ion temperature-gradient-driven modes and anomalous ion transport in tokamaks[END_REF].

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 In tokamaks, charged particles are confined by a combination of toroidal and poloidal magnetic fields that generate nested magnetic surfaces of constant magnetic flux needed for plasma confinement. Passing particles (blue trajectory) circulate along magnetic field lines in the toroidal direction ϕ, while trapped particles (red trajectory) periodically reverse their parallel velocity. For simplicity, the rapid gyromotion of charged particles is omitted and only the guiding center motion is shown. In this paper, the local transport analysis focuses at 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 inset at the top on the left 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 figure at the top on the right shows a poloidal cross-section of the plasma in the (R,Z) plane, where R is the radial distance from the torus centre and Z is the vertical coordinate; the flux surface for ρtor = 0.23 (green) and ρtor = 1 (black) are also illustrated.

field of research showed that ion scale turbulence can be partly reduced in the presence of fast ions with energies of ∼100 keV [START_REF] Citrin | Nonlinear stabilization of tokamak microturbulence by fast ions[END_REF][START_REF] Garcia | Key impact of finite-beta and fast ions in core and edge tokamak regions for the transition to advanced scenarios[END_REF][START_REF] Di Siena | Fast-ion stabilization of tokamak plasma turbulence[END_REF][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.

Simultaneously with plasma heating, alpha particles can also excite a range of instabilities in a reactor. Fast-ion-driven instabilities are often considered detrimental for plasma confinement. Not only can they 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], some of the fast-ion-driven modes can lead to enhanced 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 ion thermal improved 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 by the nonlinear interaction between MeV-range alphas, AEs destabilized by these alphas and microturbulence, studying this complex interplay at ITER-relevant conditions -both experimentally and theoretically -becomes an urgent necessity for future progress in fusion science.

For the first time, we identify in this paper reproducible experimental conditions to suppress ion scale turbulence in fusion-reactor-grade plasmas, confirmed by modeling and supported by plasma fluctuation measurements. Unlike previous fast-ion-turbulence studies [START_REF] Citrin | Nonlinear stabilization of tokamak microturbulence by fast ions[END_REF][START_REF] Garcia | Key impact of finite-beta and fast ions in core and edge tokamak regions for the transition to advanced scenarios[END_REF][START_REF] Di Siena | Fast-ion stabilization of tokamak plasma turbulence[END_REF][START_REF] Di Siena | Electromagnetic turbulence suppression by energetic particle driven modes[END_REF], turbulence suppression is now realized in JET plasmas with strong electron heating from MeV-range fast ions and simultaneously destabilized AEs. The experimental identification of novel operational conditions and the detailed theoretical understanding of the origin of turbulence suppression in these plasmas with MeV ions paves the way towards an enhanced performance of future fusion reactors with strong alpha particle heating.

Improved confinement in JET plasmas with a large population of MeV-range fast ions

The impact of fast ions on plasma dynamics was studied in JET D-3 He plasmas (n( 3 He)/n e ≈ 20-30%). Although 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], we show that the energy confinement can be actually improved in the presence of these modes. This is illustrated by comparing two JET pulses performed at the same plasma conditions (L-mode, B 0 = 3.7 T, I p = 2.5 MA, n e0 ≈ 6 × 10 19 m -3 ), differing only in the auxiliary heating scheme.

In pulse #94704 (blue lines in Figure 2), Neutral Beam Injection (NBI) was the only heating system (P aux = 12 -16 MW), injecting fast D ions with energies up to 100 keV. The measured temperature profiles of electrons and ions in the phase with P aux = 14 MW (t = 9.0 s) are shown in Figures 2(f) and (g), together with NBI power deposition profiles to bulk ions and electrons as computed by the Monte-Carlo orbit following code NUBEAM [START_REF] Pankin | The tokamak Monte Carlo fast ion module NUBEAM in the National Transport Code Collaboration library[END_REF] implemented in the modelling suite TRANSP [START_REF] Hawryluk | An empirical approach to tokamak transport[END_REF][START_REF] Ongena | Numerical Transport Codes[END_REF]. The moderately energetic D ions from NBI deposit most of their energy to bulk ions, leading to T i /T e ≈ 1.4 measured at ρ tor = 0.2. Under these conditions, the injected fast NBI ions are sub-Alfvénic (v NBI /v A ≈ 0.4) and no AEs were observed in this plasma. Red traces correspond to pulse #94701, blue traces to pulse #94704. The vertical lines show the times used in (f) and (g) for the electron and ion temperature profile comparison. The chosen timings (t = 9.5 s for #94701 and t = 9.0 s for #94704) are taken at the same total input power (14 MW) and very similar plasma densities for both discharges. In panels (f) and (g), also the deposition profiles of the total heating power on electrons and thermal ions as computed by TRANSP for both pulses are shown. In the power deposition calculations both the externally heating power (NBI + ICRF) and the electron/ion energy exchanges are considered. The error bars in (f) and (g) represent a combination of the systematic diagnostic uncertainties and time-averaged (±0.1 s) statistical variations of the measurements.

In pulse #94701 (red lines in Figure 2), D ions from NBI were accelerated to higher energies with waves in the ion cyclotron range of frequencies (ICRF). The three-ion scenario using a combination of 8 MW of NBI and 6 MW of ICRF produced deuterons with tail energies up to ∼2-3 MeV in the plasma core [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]. This is demonstrated in Figure 3(a), where measurements from the TOFOR spectrometer [START_REF] Eriksson | Measuring fast ions in fusion plasmas with neutron diagnostics at JET[END_REF] are shown. The presence of energetic D ions resulted simultaneously in long-period sawteeth and TAEs with different toroidal mode numbers (n = 2 -6), see Fig. 3(b). The observed sawtooth dynamics is similar to earlier JET experiments withICRF and monster sawteeth 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], the volume enclosed within the q = 1 surface is much smaller in these experiments in D-3 He plasmas. For example, in pulse #94701 the q = 1 surface shortly before the crash at t ≈ 9.42 s was located at ρ tor ≈ 0.2 (R inv ≈ 3.2 m), as inferred from electron cyclotron emission measurements of the T e profile. In this series of experiments, TAEs were localized at R ≈ 3.22 -3.36 m, i.e. outside the q = 1 surface, as inferred from correlation reflectometer measurements (see Fig. S2(a) in the Supplementary Information). As a result of the collisional slowing-down, highly energetic deuterons heat predominantly electrons, giving rise to a strongly peaked T e profile, see Fig. 2(f). Note that T e0 responds almost immediately to the switch-on of the ICRF power, even before TAEs are destabilized. Surprisingly, T i is also higher in pulse #94701 although the power deposition to bulk ions in the core is lower because of the reduced NBI power, compared to pulse #94704. The computed power densities shown in Figures 2(f) and (g) ignore 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 AEs (> 1 MHz), their potential impact in the reported plasmas is difficult to assess. Electron heating has been always considered a drawback with respect to ion heating by NBI as it also leads to a reduction of well-known mechanisms stabilizing ITG turbulence such as the ratio T i /T e and the E × B shearing rate (γ E×B ) [START_REF] Doyle | Chapter 2: Plasma confinement and transport[END_REF]. Yet, higher thermal confinement is obtained for pulse #94701 compared to #94704. We note that these two pulses also exhibit a large difference in the fast ion characteristics, as can be seen by comparing the measured neutron rate from D-D fusion reactions at the same amount of auxiliary heating power: 8 -9 × 10 15 s -1 in #94701 vs. 1.2 × 10 15 s -1 in #94704, see Figure 2(e). This is further evidenced in Figure 3(a) that illustrates neutron energy spectra in these two pulses at P aux = 14 MW, measured by the time-of-flight neutron spectrometer TOFOR [START_REF] Eriksson | Measuring fast ions in fusion plasmas with neutron diagnostics at JET[END_REF]. Note that shorter times-of-flight correspond to higher neutron energies, E n ≈ 2.9 MeV×(60ns/t TOF ) 2 and this reflects the presence of high-energy D ions in the plasma (cf. Figure 3(b) in [START_REF] Eriksson | Measuring fast ions in fusion plasmas with neutron diagnostics at JET[END_REF]). The observation of neutrons with t TOF = 40 -45 ns is a direct experimental confirmation for a significant amount of D ions with energies up to ∼2-3 MeV in pulse #94701. The large difference in the fast-ion content between pulses #94701 and #94704 was reproduced in numerical simulations. For pulse #94704 with ∼100 keV NBI ions, TRANSP shows T fast /T e ≈ 8.7 and n fast /n e ≈ 6%, comparable to previous studies at JET that demonstrated the beneficial effect of moderately energetic fast ions to reduce turbulence [START_REF] Citrin | Nonlinear stabilization of tokamak microturbulence by fast ions[END_REF][START_REF] Garcia | Key impact of finite-beta and fast ions in core and edge tokamak regions for the transition to advanced scenarios[END_REF][START_REF] Di Siena | Electromagnetic turbulence suppression by energetic particle driven modes[END_REF]. However, with MeV-range ions, these values become T fast /T e ≈ 33.6 and n fast /n e ≈ 3% and are much closer to those expected for fusion-born alpha particles in ITER D-T plasmas [START_REF] Garcia | Isotope and fast ions turbulence suppression effects: Consequences for high-β ITER plasmas[END_REF].

Numerical modeling of turbulence in the presence of MeV-range ions

Extensive numerical analyses of turbulence in JET pulse #94701 with a large population of MeV-range ions were carried out with the state-of-the-art gyrokinetic code Gene [START_REF] Jenko | Electron temperature gradient driven turbulence[END_REF]. Applied in its local (also called 'flux-tube') version, the Gene code has been successfully validated and has reproduced turbulence properties for various tokamak plasma experiments (see e.g. Refs. [START_REF] Citrin | Nonlinear stabilization of tokamak microturbulence by fast ions[END_REF][START_REF] Garcia | Key impact of finite-beta and fast ions in core and edge tokamak regions for the transition to advanced scenarios[END_REF][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]). In this paper, as the simulation domain a flux-tube around ρ tor = 0.23 (R = 3.25 m) was chosen, as this region is characterized by large temperature gradients, significant population of fast ions and TAEs. Experimental data at t = 9.5 s were taken as input for Gene simulations (see Table 1 in the Methods section). 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 ionsat 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 AEs 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 3(b)), these high-frequency modes appearing in the Gene simulations are identified as fast-ion-driven TAEs.

Suppression of the electrostatic turbulence in the presence of fast-ion-driven TAEs

To study the impact of MeV-range ions and fast-ion driven modes on the turbulence dynamics, nonlinear simulations with Gene were performed. Figure 4(a) shows the energy diffusivity attributed to electrostatic (ES) fluctuations for both thermal ion species, χ ES (D) and χ ES ( 3 He). The simulations reveal a strong stabilizing effect of the fast ions on the diffusivity of the thermal ions e.g. at R/L p FD = 10.7 and R/L p FD = 16.2, when the fast-ion pressure gradient is exceeded such that fast-ion driven TAEs are excited in the nonlinear turbulence regime (see also Figures 4(f) and (g)). In the presence of these modes, the computed thermal diffusivities of both thermal ion species were reduced by more than 95% compared to the values obtained without fast ions. The strongly reduced values of χ ES predicted by Gene are in a very good agreement with the estimated diffusivity of thermal ions, as calculated by the power balance analysis in TRANSP based on the experimentally measured temperature profiles (black dashed line in Figure 4(a)).

Figures 4(b)-(d) show the amplitudes of the Fourier components of the perturbed electrostatic potential φ (averaged over the radial and parallel directions) as a function of frequency and k y ρ s for three values of R/L p FD . The frequency dependence of φ additionally averaged over k y ρ s is illustrated in Figures 4(e)-(g). The turbulence spectra reveal that for R/L p F D = 5.9, the dominant modes are characterized by low frequencies and peak around k y ρ s ≈ 0.2 (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 the turbulence structure is substantially modified and the dominant modes are shifted towards higher frequencies f ≈ 200 kHz and even lower k y ρ s ≈ 0.025 -0.05. As discussed above, the characteristics of ∼200 kHz modes appearing in nonlinear Gene simulations match closely those of the experimentally observed TAEs.

Furthermore, modes at f ≈ 20 -30 kHz were observed experimentally in JET pulse #94701. Modes at similar frequencies also appeared in the Gene simulations at R/L p F D = 16.2. The bispectral analysis on the electrostatic potential computed by Gene at this value of R/L p F D 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, at a lower value R/L p F D = 10.7 such a mode coupling is not observed, while the turbulence suppression is already present in Gene simulations. Therefore, we exclude f ≈ 20 -30 kHz modes as a main cause for the observed confinement improvement.

With the knowledge of the Fourier components, the perturbed electrostatic potential φ can be calculated in real space. The surprising results, shown in Figures 5(a) and (b), illustrate the dependence of φ on the radial (x/ρ s ) and binormal (y/ρ s ) coordinates. In the absence of fast ions, one can see small-scale turbulent eddies predominately elongated along the radial direction, typical for ITG-induced transport. In contrast, intense poloidally oriented and radially sheared zonal flows appear in simulations with R/L p F D = 16.2. In this case, the zonal flow shearing rate reaches γ E×B, zonal ≈ 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. As shown in Figure 5(d), the dominant zonal structures along the y-direction are also observed for the perturbed magnetic potential A || 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]. Under these conditions and contrary to previous JET studies with moderately energetic fast ions [START_REF] Citrin | Electromagnetic stabilization of tokamak microturbulence in a high-β regime[END_REF], the electron magnetic flux remains low (χ EM (el.) ≈ 3.5 m 2 /s), despite the large-amplitude fluctuations of A induced by the fully destabilized fastion-driven modes. This is consistent with the fact that T e0 significantly increases during JET pulse #94701 unlike, e.g., in NSTX experiments [START_REF] Gorelenkov | Anomalous electron transport due to multiple high frequency beam ion driven Alfvén eigenmodes[END_REF], where the electron thermal confinement was severely degraded in the presence of AEs, thereby leading to the clamping of T e with increasing input heating power (see more details in the dedicated section of Supplementary Information).

Experimental evidence for minimized density fluctuations

For the first time, 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 reported for both electrostatic and electromagnetic perturbations in realistic experimental conditions. This coupling is further highlighted by additional bispectral analysis, suggesting a net energy transfer from TAEs to the zonal modes (see Supplementary Figure S7). The robustness of Gene results was checked against the uncertainty in the input ion temperature gradient, arising from a lack of T i data for ρ tor < 0.2 (see also Methods and Supplementary Information). Despite all the unavoidable assumptions and simplifications in the adopted local approach (e.g., for fully accurately describing zonal structures and low-n AEs radially), the mechanism is reproducible and consistent with a range of experimental observations. Modeling predictions of the turbulence suppression in these JET plasmas are supported by additional experimental evidence. Supplementary Figures S2(b) and(c) compare the level of density fluctuations measured by the correlation reflectometer in JET pulse #95669, performed at very similar operational conditions as pulse #94701. Much smaller fluctuation amplitudes are observed when MeV-range fast ions and fast-ion-driven TAEs are present in the plasma, as compared to the NBI-only phase with 100 keV fast ions. Such a reduction is also observed at R = 3.3 m, i.e. well outside the q = 1 surface, further supporting that sawtooth stabilization itself is not the main cause of the observed improved confinement.

As already outlined in the introduction, 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 [START_REF] Citrin | Nonlinear stabilization of tokamak microturbulence by fast ions[END_REF][START_REF] Garcia | Key impact of finite-beta and fast ions in core and edge tokamak regions for the transition to advanced scenarios[END_REF][START_REF] Di Siena | Fast-ion stabilization of tokamak plasma turbulence[END_REF][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. 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, the strongest reduction in the density fluctuations was measured in the same JET pulse when MeV-range ions and fast-ion-driven TAEs were present, obtained by increasing ICRF power to 6 MW. Supplementary Figure S3(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 novel 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.

Methods

GENE simulation model and input parameters

In its flux-tube version, the Gene code solves the nonlinear gyrokinetic Vlasov equations [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF] 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 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. [START_REF] Jenko | Electron temperature gradient driven turbulence[END_REF][START_REF] Merz | Gyrokinetic simulation of multimode plasma turbulence[END_REF].

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. [START_REF] Merz | Gyrokinetic simulation of multimode plasma turbulence[END_REF], which are available on the Gene website. The input parameters that have been employed in the numerical analyses are reported in Table 1.

3.00

Here, represents the inverse aspect ratio, n the density, R/Ln,T the normalized logarithmic density and temperature gradient, βe the electron-beta, and ν * the normalized collision frequency. We also report the on-axis magnetic field strength B0, the local (at ρ = 0.23) electron temperature Te and density ne, and the major radius R0. 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 4.

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., [START_REF] Görler | Multiscale effects in plasma microturbulence[END_REF].

Interpretive integrated modelling framework

Both pulses #94701 and #94704 were analysed through interpretive simulations performed with the TRANSP modelling suite [START_REF] Ongena | Numerical Transport Codes[END_REF] coupled with external heating modules NUBEAM (NBI) [START_REF] Pankin | The tokamak Monte Carlo fast ion module NUBEAM in the National Transport Code Collaboration library[END_REF] and TORIC (ICRF) [START_REF] Brambilla | Numerical simulation of ion cyclotron waves in tokamak plasmas[END_REF], and prepared with the OMFIT integrated modelling platform [START_REF] Grierson | Orchestrating TRANSP simulations for interpretative and predictive tokamak modeling with OMFIT[END_REF]. 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, and 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 Figure 1: Schematic view of the tokamak geometry.In tokamaks, charged particles are confined by a combination of toroidal and poloidal magnetic fields that generate nested magnetic surfaces of constant magnetic flux needed for plasma confinement. Passing particles (blue trajectory) circulate along magnetic field lines in the toroidal direction ϕ, while trapped particles (red trajectory) periodically reverse their parallel velocity. For simplicity, the rapid gyromotion of charged particles is omitted and only the guiding center motion is shown. In this paper, the local transport analysis focuses at 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 inset at the top on the left 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 figure at the top on the right shows a poloidal cross-section of the plasma in the (R,Z) plane, where R is the radial distance from the torus centre and Z is the vertical coordinate; the flux surface for ρtor = 0.23 (green) and ρtor = 1 (black) are also illustrated.

Figure 2 :

 2 Figure 2: Overview of the main parameters in JET pulses #94701 and #94704. The panels show (a) P NBI and P ICRF together with the total auxiliary heating power, Paux = P NBI + P ICRF ; (b) core line integrated electron density ( nedl), (c) core electron temperature Te0, (d) plasma stored energy Wp and (e) neutron rate from D-D fusion reactions.Red traces correspond to pulse #94701, blue traces to pulse #94704. The vertical lines show the times used in (f) and (g) for the electron and ion temperature profile comparison. The chosen timings (t = 9.5 s for #94701 and t = 9.0 s for #94704) are taken at the same total input power (14 MW) and very similar plasma densities for both discharges. In panels (f) and (g), also the deposition profiles of the total heating power on electrons and thermal ions as computed by TRANSP for both pulses are shown. In the power deposition calculations both the externally heating power (NBI + ICRF) and the electron/ion energy exchanges are considered. The error bars in (f) and (g) represent a combination of the systematic diagnostic uncertainties and time-averaged (±0.1 s) statistical variations of the measurements.

Figure 3 :

 3 Figure 3: Evidence of MeV-range fast D ions and fast-ion driven TAEs. Panel (a) shows TOFOR [32] measurements of the time-of-flight of the neutrons for both #94701 and #94704, indicating the presence of MeV-range fast D ions generated by the three-ion scheme. The time-of-flight t TOF of the D-D fusion-born neutrons can be related to the energy of the reactant fast D ions by kinematics considerations [36]. The error bars represent the square root of the number of counts in each channel and arise from the underlying Poisson statistics of the neutron detection process. In (b), the Mirnov coil spectrograms for pulse #94701. Strong TAE activity in the range 185-235 kHz during the application of the three-ion heating scheme is observed. Other Alfvén eigenmodes, such as EAEs and RSAEs, have been detected as well [30]. Toroidal mode numbers are indicated next to the respective traces of the detected magnetic fluctuations in the spectrogram.

Figure 4 :

 4 Figure 4: Thermal ion transport suppression simultaneous with the destabilization of FI-driven modes. In panel (a), the electrostatic thermal diffusivities of D and 3 He as a function of the fast ion normalized logarithmic pressure gradient R/Lp FD are shown. 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 of the scan can be attributed to the high-frequency and large-amplitude modulations of the energy fluxes due to strong FI-driven fluctuations of the electrostatic potential. The D and 3 He thermal diffusivities for the case without fast ions are plotted in horizontal dotted lines within shaded areas representing the standard deviations. The horizontal black dashed line represents the thermal ion energy diffusivity (χ ES ≈ 0.8 m 2 /s) obtained with TRANSP, with the grey shaded area representing the averaged uncertainty in TRANSP power balances computation (as further explained in the Methods section). Note the good agreement between the two curves only when the TAEs are nonlinearly destabilized. Panels (b-d) show the frequencies of the Fourier transform of the gyroaveraged perturbed electrostatic potential fluctuation ( φ) as a function of the binormal wavenumber and the frequency and averaged over the other spatial directions (x and z). The three cases refer to the points of the scan displayed in panel (a). In (e-g), φ is further averaged over the binormal direction. The grey shaded areas in panels (f) and (g) represent the experimentally measured range of TAE frequencies showing good agreement with the frequency peaks.

Figure 5 :

 5 Figure 5: Evidence for the nonlinear generation of zonal structures. Contour plots of the gyroaveraged perturbed electrostatic potential φ and the component in the direction of the magnetic field of the potential vector A || in the (x, y) space for (a) and (c) the case without fast ions and (b) and (d) the R/Lp F D = 16.2 case. A time average is performed for both field fluctuations at the the low field side mid-plane. The φ contour plots (a) and (b) are both normalized to the maximum (absolute) value of the R/Lp F D = 16.2 case. The same procedure is applied for A || in plots (c) and (d). This procedure highlights the increase of the intensity of both field perturbations when fast-ion driven TAEs are unstable (for R/Lp F D = 16.2) with respect to the simulation without fast ions. In (e), the electrostatic (ES) and electromagnetic (EM) energy flux spectra are shown for both electrons and 3 He with R/Lp F D = 16.2. Although the ky simulation domain extends up to kyρs = 1.2, the plot is resized for displaying properly the interesting region of the spectra.

Table 1 :

 1 Employed plasma parameters in Gene simulations modelling JET pulse #94701 at ρtor = 0.23 and t = 9.5 s. He n FD /n e T FD /T e β e[%]

		q	ŝ	T i /T e	R/L ne	R/L Te,i	n D /n e	n3 He /n e	R/L n D
	0.31	1.1	0.63	1.0	4.50	10.30	0.43	0.27	3.70
	R/L n3 ν *	B 0 [T] T e [keV]	n e [m -3 ]	R 0 [m]
	4.97	0.03	33.8	0.68	9.4 × 10 -5	3.68	4.41	5.18 × 10 19

Acknowledgements

We would like to thank Tobias Görler for providing an essential advice to ensure the correct numerical setup for the gyrokinetic simulations reported in this paper. S.M. would also like to thank Elena de la Luna for priceless support in detailing the experimental diagnostics of JET, Yann Camenen, Xavier Garbet and Andreas Bierwage for fruitful discussions about the gyrokinetic analyses, Gerardo Giruzzi for valuable suggestions on the article strategy, and Aaron Ho for assisting the author in processing the experimental data. We sincerely thank our colleagues Matteo Baruzzo and Filomena Nave for the excellent preparation and execution of the experiments discussed in this paper. 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.

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

S.M. performed the gyrokinetic simulations and the subsequent analyses, with the fundamental support of J.G., D.Z. and S.B. The analysis shedding lights on the excitation of the fast-ion driven modes were performed by D.Z. The reported experiments were devised and jointly led by Ye.O.K., J.O., J.G. and M.N., with the key assistance of M.D. in analyzing the experimental outcomes. Crucial input data for the gyrokinetic analyses were provided by Ž. Š., G.S., Ye.O.K. and M.D., whereas J.E. and A.S. supported the analysis of the integrated modelling results. The manuscript was written by S.M., J.G., D.Z, Ye.O.K. and J.O. with feedback by all the authors.