Francis Filbet

Charles Prouveur

HIGH ORDER TIME DISCRETIZATION FOR BACKWARD SEMI-LAGRANGIAN METHODS

Keywords: High order time discretization, Semi-Lagrangian scheme, Vlasov-Poisson model, Guidingcentre model, Plasma physics

We introduce different high order time discretization schemes for backward semi-Lagrangian methods. These schemes are based on multi-step schemes like Adams-Moulton and Adams-Bashforth schemes combined with backward finite difference schemes. We apply these methods to transport equations for plasma physics applications and for the numerical simulation of instabilities in fluid mechanics. In the context of backward semi-Lagrangian methods, this time discretization strategy is particularly efficient and accurate when the spatial error discretization becomes negligeable and allows to use large time steps.

Introduction

A model which can be used in many cases for the study of plasma as well as beam propagation is the Vlasov equation coupled with the Maxwell or Poisson equations to compute the self-consistent fields. It describes the evolution of a system of particles under the effects of external and self-consistent fields. The unknown f (t, x, v), depending on the time t, the position x, and the velocity v, represents the distribution of particles in phase space for each species with (x, v) ∈ R d × R d , d = 1, .., 3. Its behaviour is given by the Vlasov equation, (1.1) ∂f ∂t

+ v • ∇ x f + F(t, x, v) • ∇ v f = 0,
where the force field F (t, x, v) is coupled with the distribution function f giving a nonlinear system. We first define ρ(t, x) the charge density and J(t, x) the current density which are given by

ρ(t, x) = q R d f (t, x, v)dv, J(t, x) = q R d v f (t, x, v)dv,
where q is the elementary charge. For the Vlasov-Poisson model

(1.2) F(t, x, v) = q m E(t, x), E(t, x) = -∇ x φ(t, x), -∆ x φ = ρ ε 0 ,
where m represents the mass of one particle. On the other hand for the Vlasov-Maxwell model, we have

F(t, x, v) = q m (E(t, x) + v × B(t, x) ),
and E, B are solutions of the Maxwell equations

                   ∂E ∂t -c 2 ∇ × B = - J ε 0 , ∂B ∂t + ∇ × E = 0, ∇ • E = ρ ε 0 , ∇ • B = 0,
with the compatibility condition ∂ρ ∂t + div x J = 0, which is verified by the Vlasov equation solution.

Later on we will also consider the guiding centre model, which has been derived to describe the evolution of the charge density in a highly magnetized plasma in the transverse plane of a tokamak. This model is described as follows: we now consider the density ρ solution to (1.3)

   ∂ρ ∂t + U • ∇ρ = 0, -∆φ = ρ,
where the velocity U = (-∂ y φ, ∂ x φ) is divergence free. The numerical resolution of the Vlasov equation and related models is usually performed by Particle-In-Cell (PIC) methods which approximate the plasma by a finite number of particles. Trajectories of these particles are computed from characteristic curves given by the Vlasov equation, whereas self-consistent fields are computed on a mesh of the physical space. This method yields satisfying results with a relatively small number of particles but it is subject to fluctuations, due to the numerical noise, which are difficult to control.

To develop accurate and stable numerical techniques for plasma turbulence, numerical methods based on direct numerical simulation techniques have been developped. The Vlasov equation is discretized in phase space using either semi-Lagrangian [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF][START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: theoretical analysis and application to the Vlasov-Poisson system[END_REF][START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF], finite element [START_REF] Guo | Hybrid semi-Lagrangian finite element-finite difference methods for the Vlasov equation[END_REF], finite difference [START_REF] Umeda | A non-oscillatory and conservative semi-Lagrangian scheme with fourth-degree polynomial interpolation for solving the Vlasov equation[END_REF] or discontinuous Galerkin [START_REF] De Dios | Discontinuous Galerkin methods for the multidimensional Vlasov-Poisson problem[END_REF][START_REF] Heath | A discontinuous Galerkin method for the Vlasov-Poisson system[END_REF] schemes.

Among them, the semi-Lagrangian method consists of computing the distribution function at each grid point by following the characteristic curve ending there. Then to compute the value of the distribution function at the origin of the characteristic, a high order interpolation method is needed.

Concerning the time discretization of the backward differential system corresponding to the characteristic curved, it is often performed via a second order splitting scheme, but higher order techniques are also available [START_REF] Afeyan | Simulations of Kinetic Electrostatic Electron Nonlinear (KEEN) Waves with Two-Grid, Variable Velocity Resolution and High-Order Time-Splitting, EPJD[END_REF]. This strategy is particulary interesting when each step of the splitting can be explicitly solved. However, it is not always possible to apply a splitting strategy, hence the time discretization technique is often limited to a leap-frog scheme or a second order predictor-corrector scheme. The aim of the present work is to propose and validate high order techniques to preserve the order of accuracy of the spatial discretization.

This paper is organized as follows : in Section 2 we recall the backward semi-Lagrangian method then in Section 3 we develop different time discretization schemes, based on multi-step methods. Both classical second order schemes and "new" third and fourth order schemes will be proposed. In Section 4, we present several numerical results. In the first, a two-dimensional transport equation with gaussian initial data is investigated to recover the expected order of accuracy. In the second, we perform numerical simulations on the simplified paraxial Vlasov-Poisson model in two dimensions. The last three tests are applied to the guiding centre model to compare the performances of the schemes and to investigate the ability to recover stability results and the development of instabilities.

The backward semi-Lagrangian method

To construct a semi-Lagrangian method for (1.1) or (1.3), we consider a generic transport equation written in an advective form

(2.1) ∂f ∂t + A • ∇f = 0,
where

A : R + × R d → R d , d ≥ 1
and a set of mesh points (x i ) i∈Z Z .The main feature of the backward semi-Lagrangian method is to compute an approximated solution to (2.1) by following the characteristics backwards in time

(2.2)      dX dt = A, X(t n+1 ) = x i .
where x i is a mesh point. Thus, we evaluate the distribution function f at this point at time t n+1 by following its characteristic curve backward in time until we reach the time when the values of the distribution function f are known, t n or t n-1 . In general the coordinates of the characteristic at this time do not coincide with the mesh points, therefore we need to interpolate the value of the distribution function at this point. Hence, this method can be divided in two steps :

• the first is the resolution of the characteristic differential equation (2.2). This system's solution at time t n , respectively t n-1 depending on the time scheme used, it will be denoted by X n , respectively X n-1 ; • the second step is the interpolation in order to compute an approximation of f (t n , X(t n )), since the value of the distribution function is constant along the characteristic, we denote by f n+1 the approximation of f (t n+1 ) given by

f n+1 (x i ) = Π h f n (X n ), i ∈ ZZ,
where Π h is an interpolation operator. A special case of semi-Lagrangian method for Vlasov-Poisson system, first introduced by Cheng and Knorr [START_REF] Shoucri | Numerical integration of the Vlasov equation[END_REF], is based on a time splitting which enables an exact computation of the characteristics coupled with a cubic spline interpolation for the phase space reconstruction. This algorithm was subsequently applied in many plasma physics papers, see for example [START_REF] Feix | Eulerian codes for the Vlasov equation Series[END_REF][START_REF] Ghizzo | A Vlasov code for the numerical simulation of stimulated Raman scattering[END_REF] and references therein. This method was cast into the more general framework of semi-Lagrangian methods by E. Sonnendrücker et al. [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF], and successfully adapted to beam physics problems, namely the simulation of space charge waves. For Vlasov type equations, one of the most popular interpolation is the multi-dimensional cubic spline interpolation which is obtained thanks to the tensor product of cubic splines [START_REF] Shikin | Handbook on splines for the user[END_REF]. Another flavor of the semi-Lagrangian method was introduced by Nakamura and Yabe and called the Cubic Interpolated Propagation (CIP) method. It is based on a Hermite interpolation for which the gradients of the distribution function are also advanced along the characteristics [START_REF] Nakamura | Cubic interpolated propagation scheme for solving the hyper-dimensional Vlasov-Poisson equation in phase space[END_REF]. It needs the storage of f , ∇ x f , and ∇ v f , therefore in order to reduce memory consumption, the mesh on which it is applied should be coarser. We refer to [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations[END_REF] for a comparison of different reconstruction techniques.

Most of these methods based on a time splitting discretization are particularly efficient for classical Vlasov-Poisson or Vlasov-Maxwell systems since the characteristic curves corresponding to the split operator simply become straight lines and can be solved exactly. Therefore, the numerical error is only due to the splitting in time and the phase space discretization of the distribution function. Furthermore for such time splitting schemes, the semi-Lagrangian methods on Cartesian grids coupled with Lagrange, Hermite or cubic spline interpolation techniques are conservative [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF].

However, for more elaborated kinetic equations, or even the two dimensional guiding centre model [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF], time splitting techniques cannot necessarily be applied. Thus characteristic curves are more sophisticated and require a specific time discretization. For instance, in [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF] several numerical solvers have been developed using an Eulerian formulation for gyro-kinetic models. Assume that for all time t m , m ≤ n, approximated solutions of the distribution f m and of the velocity field A m are known. The semi-Lagrangian technique based on the leap-frog scheme (2.2) consists in computing X n-1 such that

       X n+1 -X n-1 2∆t = A n (X n ), X n+1 = x i ,
where the two unknowns X n-1 and X n are related by the mid-point rule

X n+1 + X n-1 2 = X n .
Finally, the unknown X n-1 is solution to the following nonlinear problem:

(2.3)      X n+1 -X n-1 2∆t = A n X n+1 + X n-1 2 , X n+1 = x i .
However, this scheme is known to generate spurious oscillations in the non-linear phase when small structures occur and it is difficult to distinguish physical and numerical oscillations. Moreover, for these models semi-Lagrangian methods are no more conservative, hence the long time behavior of the numerical solution may become unsuitable.

An alternative time discretization consists in applying a predictor-corrector technique [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF]. Assume that at time t n , approximated solutions of the distribution f n and of the velocity field A n are known. The prediction stage is performed with a half time step : a first order approximation is applied and X n is given by

(2.4)      2 X n+ 1 2 -X n ∆t = A n (X n ), X n+ 1 2 = x i .
Then we compute the value of the distribution function at time

t n+ 1 2 = t n + ∆t 2 for each mesh point x i , i ∈ ZZ, (2.5) f n+ 1 2 (x i ) = Π h f n (X n ).
For nonlinear problems such as (1.1) or (1.3) the velocity field A depends on f , hence its approximation A n+1/2 at time t n+ 1 2 can be computed from f n+1/2 . The third stage consists in computing a new approximation of X n thanks to the previous prediction and evaluation at time t n+1/2 with a second order leap-frog scheme

(2.6)      X n+1 -X n ∆t = A n+ 1 2 X n+1 + X n 2 ,
X n+ 1 2 = x i . Finally, for the fourth stage, the value of the solution is given by a last interpolation :

(2.7)

f n+1 (x i ) = Π h f n (X n ).
This method gives much more stable and accurate results than the previous leap-frog scheme (2.3), but it requires one more interpolation and thus has a higher computational cost compared to other methods of the same order of precision. Furthermore, it is limited to second order accuracy whereas the interpolation procedure commonly based on cubic spline techniques of high degree interpolations is third order accurate.

For this purpose, we want to develop a class of high order numerical schemes specifically designed for backward semi-Lagrangian method. In the following we present several schemes based on Adams-Moulton and Adams-Bashforth time discretization techniques applied in the framework of backward semi-Lagrangian methods.

High order time discretization for Backward Semi-Lagrangian methods

In this section, we fix a time step ∆t > 0 and assume that the velocity field A m and the distribution function f m have been approximated for all time step t m = m ∆t, m ≤ n, we look for an approximation X n of the solution the system (2.2) at time t n .

Hence, we will construct different high order scheme in the following form (3.8)

X n = x i -∆t Φ(∆t, t n+1 , x i ),
where Φ also depends on A at different times. This framework allows to study consistency and stability of the different numerical schemes using a systematic way. Indeed, for the consistency analysis we have [START_REF] Hairer | Solving ordinary differential equations. I. Nonstiff problems[END_REF] Proposition 3.1. Assume that the velocity field A is given. The scheme (3.8) is consistent if we have Φ(0, t n+1 , x i ) = A(t n+1 , x i ). Moreover, the scheme is consistent of order p if

∂ k Φ ∂s k (s = 0, t n+1 , x i ) = (-1) k+1 k + 1 A (k) (t n+1 , x i ), 0 ≤ k ≤ p,
where A (0) = A and for k ≥ 1

A (k) = ∂ ∂t A (k-1) + A • ∇ x A (k-1) .
3.1. Second order Adams-Moulton scheme. We first consider the implicit Adams-Moulton scheme, which is second order accurate :

X n+1 -X n ∆t = 1 2 (A n+1 (X n+1 ) + A n (X n )).
Unfortunately, the velocity field A is not known at time (t n+1 , thus a high order reconstruction is performed using an affine extrapolation in time using previous time steps t n-1 and t n , it yields

A n+1 (X n+1 ) = 2 A n (X n+1 ) -A n-1 (X n+1 ).
Therefore the approximation X n is found by solving:

(3.9)      X n+1 -X n ∆t = 1 2 2 A n (X n+1 ) -A n-1 (X n+1 ) + A n (X n ) , X n+1 = x i .
This nonlinear system can be either solved using a classical fixed point or Newton method or approximated by a linear system using a piecewise affine approximation of A n . It can be written in the form

X n = x i -∆t Φ(∆t, t n+1 , x i ),
where the function Φ is implicitly defined.

Proposition 3.2. Assume that the velocity field

A ∈ C 1 (R + × R d
) is given. The scheme (3.9) is consistent and second order accurate.

Proof. Assume that A is given and A ∈ C 1 . From the definition of the second order Adams-Moulton scheme (3.9), we define implicitly Φ such that

Φ(∆t, t n+1 , x i ) = 1 2 2A(t n+1 -∆t, x i ) -A(t n+1 -2∆t, x i ) +A t n+1 -∆t, x i -∆t Φ(t n+1 , x i , ∆t) .
Therefore, we easily verify that Φ(0, t n+1 , x i ) = A(t n+1 , x i ) anf furthermore by differentiating with respect to ∆t, we get that

∂Φ ∂s (s = 0, t n+1 , x i ) = - 1 2 ∂A ∂t + A • ∇ x A (t n+1 , x i ).
Hence the scheme is second order accurate.

3.2.

Third order Adams-Moulton scheme. Now we start from the third order accurate Adams-Moulton scheme applied to (2.2), it gives

X n+1 -X n ∆t = 1 12 5 A n+1 (X n+1 ) + 8 A n (X n ) -A n-1 (X n-1 ) .
In order to define rigourously the numerical approximation X n we need to fix two different issues. The first one is to obtain an approximation A n+1 of the velocity field A at time t n+1 , which preserves the order of accuracy up to third order. Therefore, as for the previous scheme, we extrapolate the velocity field A from the values at the previous time steps using a Lagrange interpolation polynomial of degree two, that is,

A n+1 (X n+1 ) = 3 A n (X n+1 ) -3 A n-1 (X n+1 ) + A n-2 (X n+1 ).
The second issue we have to address, comes from the fact that the third order Adams-Moulton scheme is a two steps scheme i.e. we have introduced an additional unknown X n-1 . Hence, another relation between the unknowns X n-1 and X n is needed and, as for the previous extrapolation, it has to preserve the order of accuracy of the scheme. Let us consider the following polynomial P 2 (t) such that

             P 2 (t n+1 ) = X n+1 , P 2 (t n ) = X n , dP 2 dt (t n+1 ) = A n+1 (X n+1 ),
where the velocity field A n+1 has been computed using the previous extrapolation step. Finally, we define X n-1 = P 2 (t n-1 ), which yields the following relation

X n-1 = 4X n -3 X n+1 + 2 ∆t A n+1 .
Gathering the previous approximations, it gives the following one step scheme (3.10)

                     X n+1 -X n ∆t = 1 12 5 A n+1 (X n+1 ) + 8A n (X n ) -A n-1 (X n-1 ) , A n+1 (X n+1 ) = 3 A n (X n+1 ) -3 A n-1 (X n+1 ) + A n-2 (X n+1 ), X n-1 = 4 X n -3 X n+1 + 2 ∆tA n+1 (X n+1
),

X n+1 = x i .
It can be written in the form

X n = x i -∆t Φ(∆t, t n+1 , x i ),
where the function Φ is implicitly defined. By following the lines of Proposition 3.2 and after some calculus we show that Proposition 3.3. Assume that the velocity field A ∈ C 2 (R + × R d ) is given. The scheme (3.10) is consistent and third order accurate.

Proof. Assume that A is given and A ∈ C 1 . We first construct an approximation of A at time t n+1 as a function of ∆t, that is,

A n+1 (∆t) = 3 A(t n+1 -∆t) -3 A(t n+1 -2∆t) + A(t n+1 -3∆t).
From the definition of the third order Adams-Moulton scheme (3.10), we define implicitly Φ such that

Φ(∆t, t n+1 , x i ) = 1 12 5A n+1 (∆t, x i ) + 8 A t n+1 -∆t, x i -∆t Φ(∆t, t n+1 , x i ) -A t n+1 -2∆t, x i -4∆t Φ(∆t, t n+1 , x i ) + 2 ∆t A n+1 (∆t, x i ) .
Therefore, we easily verify that Φ(0, t n+1 , x i ) = A(t n+1 , x i ) anf furthermore by differentiating with respect to ∆t, we get that

∂Φ ∂s (s = 0, t n+1 , x i ) = ∂A ∂t + A • ∇ x A (t n+1 , x i ),
which proves that the scheme is second order accurate. Finally, to show that (3.10) is third order, we compute

∂ 2 Φ ∂s 2 (s = 0, t n+1 , x i ) = 1 3
∂ 2 A ∂s 2 + 2 ∂ ∂s ∇ x A + A (1) ∇ x A + A • ∇ 2 x A • A , = 1 3 A 2 (t n+1 , x i ).
3.3. Third order Adams-Bashforth scheme. We now start from the third order explicit Adams scheme applied to 2.2, it first gives that

X n+1 -X n ∆t = 1 12 23 A n (X n ) -16 A n-1 (X n-1 ) + 5 A n-2 (X n-2 ) .
The stencil here is different from the previous one for the Adams-Moulton scheme since we do not need to compute an approximation for A n+1 , we use a slightly different polynomial to approximate X n-1 and X n-2 .

The main issue comes now from the fact that the third order Adams scheme is a three steps scheme i.e. we have introduced additional unknowns X n-2 and X n-1 . To compute them, let us consider the following polynomial P 2 (t) such that

             P 2 (t n+1 ) = X n+1 , P 2 (t n ) = X n , dP 2 dt (t n ) = A n (X n ).
We set X n = x i -c n and from the polynomial P 2 , we can define X n-1 and X n-2 by

X n-1 = x i -2 ∆t A n (X n ), X n-2 = x i + 3 c n -6 ∆t A n (X n ).
Gathering these approximations, it leads to the resolution of a nonlinear problem which is again solved with a Newton algorithm:

(3.11)                              X n+1 -X n ∆t = 1 24 23 A n (X n ) -16 A n-1 (X n-1 ) + 5 A n-2 (X n-2 ) , X n-2 = x i + 3 c n -6 ∆t A n (X n ), X n-1 = x i -2 ∆t A n (X n ), X n = x i -c n , X n+1 = x i .
Let us note that this scheme requires more memory since it needs the storage of the last three time steps of the velocity field. It is also slightly more complicated than the implicit scheme of the same order.

By following the lines of Proposition 3.3 and after some calculus we show that Proposition 3.4. Assume that the velocity field A ∈ C 2 (R + × R d ) is given. The scheme (3.11) is consistent and third order accurate.

3.4.

Fourth order Adams-Moulton scheme. We first consider the fourth order Adams-Moulton scheme:

X n+1 -X n ∆t = 1 24 9 A n+1 (X n+1 ) + 19 A n (X n ) -5 A n-1 (X n-1 ) + A n-2 (X n-2 ) .
The stencil is now quite large and we need an approximation of A n+1 of order O(∆t 4 ). Therefore we take:

A n+1 (x) = 4 A n (x) -6 A n-1 (x) + 4 A n-2 (x) -A n-3 (x).
To complete the numerical approximation, new relations are needed for X n-1 and X n-2 . To this aim, we consider the following cubic Hermite interpolation polynomial P 3 (t) such that

                     P 3 (t n+1 ) = X n+1 , P 3 (t n ) = X n+1 , dP 3 dt (t n ) = A n (X n ), dP 3 dt (t n+1 ) = A n+1 (X n+1 ),
where the velocity field A n+1 has been computed using the previous extrapolation step.

Thus we set X n = x i -c n and X n+1 = x i , and obtain

X n-1 = x i + 4c n -2∆t 2A n (X n ) + A n+1 (X n+1 ) , X n-2 = x i + 9c n -3∆t 3A n (X n ) + A n+1 (X n+1 ) .
Gathering the previous approximations, it gives the following one step scheme (3.12)

                                                     X n+1 -X n ∆t = 1 24 9A(X n+1 ) + 19A n (X n ) -5A n-1 (X n-1 ) + A n-2 (X n-2 ) , A n+1 (X n+1 ) = 4 A n (X n+1 ) -6 A n-1 (X n+1 ) + 4 A n-2 (X n+1 ) -A n-3 (X n+1 ), X n-2 = x i + 27c n -6∆t 3A n (X n ) + 2A n+1 (X n+1 ) , X n-1 = x i + 4c n -2∆t 2A n (X n ) + A n+1 (X n+1 ) , X n = x i -c n , X n+1 = x i .
which can be solved through fixed point iteration procedure. By following the lines of Proposition 3.3 and after some heavy calculus we show that Proposition 3.5. Assume that the velocity field A ∈ C 3 (R + × R d ) is given. The scheme (3.12) is consistent and fourth order accurate.

3.5.

Fourth order Adams-Bashforth scheme. We now start from the fourth order explicit Adams scheme applied to 2.2, it first gives that

X n+1 -X n ∆t = 1 24 55 A n (X n ) -59 A n-1 (X n-1 ) + 37 A n-2 (X n-2 ) -9 A n-3 (X n-3 ) .
The stencil here is different from the previous one for the Adams-Moulton scheme since we do not need to compute an approximation for A n+1 , we use a slightly different polynomial to approximate X n-1 , X n-2 and X n-3 .

The main issue comes now from the fact that the fourth order Adams scheme is a four steps scheme i.e. we have introduced additional unknowns X n-3 , X n-2 and X n-1 . To compute them, let us consider the following polynomial P 3 (t) such that We proceed as for the implicit third order Adams scheme, and obtain the following two steps scheme

(3.14)                              X n+1 -X n-1 2∆t = 1 6 (A n+1 + 4 A n + A n-1 )., A n+1 (X n+1 ) = 4 A n (X n+1 ) -6 A n-1 (X n+1 ) + 4 A n-2 (X n+1 ) -A n-3 (X n+1 ), X n-1 = x i + 4c n -2∆t 2A n (X n ) + A n+1 (X n+1 ) , X n = x i -c n , X n+1 = x i .
This scheme is not a one step scheme and does not enter in the class presented at the beginning.

Numerical simulations

In this section we will compare the different time discretization for the backward semi-Lagrangian method using the classical cubic spline interpolation and also the fifth order Hermite interpolation. Reference solutions will be computed using a fifth order finite difference technique with a WENO reconstruction coupled with a fourth order Runge-Kutta scheme for the time discretization with a small time step [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations[END_REF]. The first numerical test concerns a simple linear advection problem for which we can compute an exact solution and verify the order of accuracy. The next simulations concern the development of instabilities for the paraxial Vlasov-Poisson approximation and the guiding centre model.

Test 1 : linear advection. We first consider a linear advection problem given by

     ∂f ∂t + V • ∇ x f = 0, x ∈ R 2 , t > 0, f (0, x) = e -x 2 ,
where V(t, x) = tan(t) x/2. For this equation, we can compute an exact solution, it yields

f (t, x) = e -x 2 cos(t) .
We define the discrete L 2 norm of the error as;

ε 2 = max 0≤n≤Nt i h 2 f n i -f (t n , x i ) 2 1/2 .
On the one hand, the computational domain is reduced to the square [-6, 6] 2 and we introduce the grid size h = 1/N x such that the grid point are given by x i = i h where i ∈ ZZ 2 . On the other hand, for the time discretization, we choose the time step ∆t = 1/N t and t n = n∆t, with n ∈ IN. Therefore, the mesh can be described only by the couple (N t , N x ). To estimate the order of accuracy of the schemes we compute a numerical approximation and refine the time step ∆t according to the space step h in such a way the CFL condition associated to the advection equation is equal to 2.

Absolute error in discrete L 2 norms at time T = 1 are shown in Table 1 for all the schemes presented in the previous section. In order to focus on the time discretization error, a large time step is used with a fine grid for the spatial discretization. As expected the order of accuracy is satisfied for all second, third and fourth order schemes.

Furthermore, in Table 2 we indicate the computational time for each discretization . We notice that the Adams-Moulton schemes are actually faster than the Adams-Bashforth schemes and much more precise. On the other hand, the second order Adams-Moulton is faster than the second order Predictor-Corrector scheme, but the numerical error is ten times larger. However, the third order Adams-Moulton scheme (3.10) is a good compromise between accuracy and efficiency.

(N t , N x ) (40,
4.2. Test 2 : paraxial Vlasov-Poisson model. We apply the numerical methods presented in previous sections to the following Vlasov-Poisson system satisfied by f (t, r, v), where

r ∈ R, v ∈ R [4, 11] (4.1)
         ∂f ∂t + v ε ∂f ∂r + E f - r ε ∂f ∂v = 0, 1 r ∂ ∂r (rE f ) = R f dv.
The electric field can be expressed explicitly as follows

E f (t, r) = 1 r r 0 s ρ(t, s) ds,
where ρ(t, r) = R f (t, r, v)dv, hence we will compute E f by a simple numerical integration. The initial condition is chosen as a Gaussian in velocity multiplied by a regularized step function in r:

(4.2) f 0 (r, v) = 4 √ 2πα χ(r) exp - v 2 2α , with
χ(r) = 1 2 erf r + 1.2 0.3 - 1 2 erf r -1.2 0.3
and α = 0.2. The Vlasov-Poisson system (4.1) conserves mass

d dt R 2 f (t, r, v)drdv = 0,
and also

L p norms for 1 ≤ p ≤ ∞ d dt ||f (t, r, v)|| L p (R 2 ) = 0.
Therefore, the evolution in time of these quantities will be observed for each scheme, especially the mass and the positivity. We will also investigate the time evolution of the kinetic energy of the Vlasov-Poisson system (4.1) :

(4.3) E(t) = R 2 v 2 2 f dr dv.
In the following we take the parameter ε = 0.7 and the computational domain is (r, v) ∈ Ω = [-4, 4] 2 . On the one hand, a reference solution is presented in Figure 2 using a grid of size 800 × 800 and a time step of ∆t = 1/400 using a fifth order space discretization for the transport equation coupled with a fourth order Runge-Rutta method for the time discretization and a fourth order method for the computation of the electric field.

On the other hand, a backward semi-Lagrangian method is applied with different type of time discretizations coupled with a fifth order Hermite interpolation. The numerical results, corresponding to the distribution function at different times, are presented in Figure 3. The grid is composed of 400 × 400 mesh points whereas the time step is ∆t = 1/100 for the different Adams-Moulton schemes (3.9) and (3.10) and for the leap-frog scheme (2.3).

We first present the time evolution of theoretically conserved quantities like mass, L 2 norm in Figure 1. Concerning the conservation of mass, the fourth order schemes gives clearly better accuracy than the second and third order schemes. However, for the time evolution of the L 2 norm, it is not possible to distinguish the results given by the different time solvers. Finally, in Figures 2 and3 we report the snapshots of the distribution function at different times t = 5, 15 and 50. High order Adams-Moulton schemes (3.9), (3.10) and (3.12), not shown here , give slighly the same resuts as the ones obtained with the finite difference method. On the other hand, the leap-frog scheme generates instabilities as it has been already observed in [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations[END_REF]. Here we consider the model in a disk domain

D = {x ∈ R 2 : x ≤ R}
and assume that the electric potential is vanishing at the boundary (4.5) φ(x) = 0, x ∈ ∂D.

Then assuming that we have zero flux at the boundary, the guiding centre model verifies the following properties :

(1) Positivity of density 0 ≤ ρ(t, x).

(2) Mass conservation d dt D ρ(t, x)dx = 0.

(

) L p norm conservation, for 1 ≤ p ≤ ∞ d dt ||ρ|| L p (D) = 0. 3
(4) Energy conservation 1 2

d dt D |∇φ| 2 (t, x)dx = d dt D ρ(t, x) φ(t, x)dx = 0.
For this test we consider the diocotron instability for an annular electron layer. The initial data is given by where ε is a small parameter, θ = atan(y/x). In the following tests, we take ε = 0.001, r -= 5, r + = 8, = 7, furthermore ∆t = 0.05 and h = ∆x = ∆y = 24/1000. We compare the results obtained with the backward semi-Lagrangian method presented in Figure 6, to the reference results obtained in Figure 5. Here, the different schemes give very similar results. The slight differences can be seen looking at Figure 4. Indeed mass, L 2 norm and energy conservation are much better with the third order Adams-Moulton scheme than with the second order. Oddly here the results with the Adams-Moulton scheme of order four are not better than those obtained with the third order. Regarding the minimum of the distribution function there are no significant differences. Hence, as it has been observed before, the third order Adams-Moulton scheme (3.10) is a good compromise here between accuracy and efficiency.

ρ 0 (x) = (1 + ε cos( θ)) exp (-4(r -6.5) 2 ), if r -≤ x ≤ r + , 0, otherwise, (a
4.4. Test 4 : Kelvin-Helmoltz instability in a rectangle. We still consider the guiding centre model but now the computational domain is a rectangle with periodic boundary conditions in the x direction and Dirichlet in the y direction [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF]. The initial distribution function is ρ 0 (x) = sin(y) (1 + 0.005 cos(k x x)), x = (x, y) ∈ (0, Lx) × (0, 2π), with k x = 2π/L x . We choose L x = 10 in such a way that the solution will generate an instability [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of Vlasov Equations[END_REF].

For this test, we choose to run the code on very fine meshes 2048×2048 points and the Poisson equation with periodic boundary conditions in x and homogeneous Dirichlet condition in y is dicretized using a classical fast Fourier transform giving high order accuracy. Concerning the time discretization we still apply the second (3.9), third (3.10) and fourth (3.12) order Adams-Moulton schemes In Figure 7, we present the time evolution of macroscopic quantities In this work, we are interested in the long time simulation of (4.6) for small initial perturbations of a stable steady state. Recently, J. Bedrossian and N. Masmoudi [START_REF] Bedrossian | Asymptotic stability for the couette flow in the 2D Euler equations[END_REF] proved that for a smooth and small enough initial perturbation of amplitude , the solution to (4.6) converges to a new equalibrium and that the velocity field U = (U x , U y ) satisfies when t goes to infinity

U y L 2 ≤ C t 2 and U x -< U x > L 2 ≤ C t ,
where < U x > denotes the local average of U x with respect to the x variable. Here we consider the equilibrium ρ eq = sin(y) and introduce a small perturbation δρ = ε cos(kx) sin(y). In Figure 8, we present the time evolution of these quantities with respect to time obtained with the third order Adams-Moulton scheme. The numerical solution follow the correct long time behavior and the convergence to equilibrium is given with the expected rate of convergence : 1/t for the quantity U x -< U x > L 2 whereas the quantity U y L 2 is of order 1/t 2 . Finally in Figure 9, we present the snapshots of the vorticity ρ -ρ eq for different times t = 20, 40 and 60 and observed that due to periodic boundary conditions oscillations in the y component appear and generate higher and higher frequencies. 3) (middle) a second order Adams-Moulton scheme (3.9) and (bottom) a third order Adams-Moulton scheme (3.10) using a grid of size 1000 × 1000.

Conclusion and perspective

In this paper we have developed and compared several time discretization techniques for the backward semi-Lagrangian method. We obverve that when the spatial mesh is fine enough and the spatial discretization error negligeable, the high order Adams-Moulton scheme brought better results than the classical schemes. In particular, the results reveal a major improvements in terms of mass and energy conservations. Hence the application of a higher order norm for a second order Adams-Moulton scheme (3.9), a third order Adams-Moulton scheme (3.10) and a fourth order Adams-Moulton scheme (3.12) using a grid of size 2048 × 2048.

schemes is justified, although in our numerical simulations four order schemes brought few improvements compared to the third order discretization. Finally, the third order Adams-Moulton scheme (3.10) is computationally interesting, robust and accurate for various numerical tests.

Figure 1 .

 1 Figure 1. Test 2 : paraxial Vlasov-Poisson model. Time evolution of conserved quantities (a) mass and (b) L 2 norm for a second order Adams-Moulton scheme (3.9), a third order Adams-Moulton scheme (3.10) and a fourth order Adams-Moulton scheme (3.12) using a grid of size 400 × 400.

Figure 2 . 4 . 3 . 2 ( 4

 24324 Figure 2. Test 2 : paraxial Vlasov-Poisson model. Reference solution computed with a high order finite difference method [24] using a grid of size 800 × 800.

Figure 3 .

 3 Figure 3. Test 2 : paraxial Vlasov-Poisson model. Numerical solution computed with (top) a leap-frog scheme (2.3) (middle) a second order Adams-Moulton scheme (3.9) and (bottom) a third order Adams-Moulton scheme (3.10) using a grid of size 400 × 400.

Figure 4 .

 4 Figure 4. Test 3 : diocotron instability. Time evolution of conserved quantities (a) mass and (b) energy (c) L 2 norm and (d) L 4 norm for a second order Adams-Moulton scheme (3.9), a third order Adams-Moulton scheme (3.10) and a fourth order Adams-Moulton scheme (3.12) using a grid of size 1000 × 1000.

Figure 5 . 4 . 5 .

 545 Figure 5. Test 3 : diocotron instability.Reference solution computed with a high order finite difference method[START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations[END_REF] using a grid of size 1000 × 1000.

Figure 6 .

 6 Figure 6. Test 3 : Diocotron instability. Numerical solution computed with (top) a leap-frog scheme (2.3) (middle) a second order Adams-Moulton scheme (3.9) and (bottom) a third order Adams-Moulton scheme (3.10) using a grid of size 1000 × 1000.

Figure 7 .

 7 Figure 7. Test 4 : Kelvin-Helmoltz instability in a rectangle. Time evolution of conserved quantities (a) mass, (b) energy, (c) L 2 norm and (d) L 4norm for a second order Adams-Moulton scheme (3.9), a third order Adams-Moulton scheme (3.10) and a fourth order Adams-Moulton scheme (3.12) using a grid of size 2048 × 2048.

Figure 8 .Figure 9 .

 89 Figure 8. Test 5 : shear flow for 2D Euler equations. Time evolution of U x -< U x > L 2 and U y L 2 using the third order Adams-Moulton scheme (3.10) with a grid of size 256 × 4096.

Table 1 .

 1 Comparison of different time discretization schemes for the backward semi-Lagrangian method : error L 2 -norm and order of convergence r.

		200)	(80,400)	(160,800)
	• 2	r	• 2	r	• 2	r

Table 2 .

 2 Comparison of computational time for the different schemes (in seconds)

We set X n = x i -c n and from the polynomial P 3 , we can define X n-1 , X n-2 and X n-3 by

Gathering these approximations, it leads to the resolution of a nonlinear problem which is again solved with a Newton algorithm:

Let us note that this scheme requires more memory since it needs the storage of the last four time steps of the velocity field. It is also slightly more complicated than the implicit scheme of the same order.

By following the lines of Proposition 3.3 and after some heavy calculus we show that Proposition 3.6. Assume that the velocity field A ∈ C 3 (R + × R d ) is given. The scheme (3.13) is consistent and fourth order accurate.

3.6. Milne-Simpson scheme. The last scheme we want to investigate is based on the two steps fourth order Milne-Simpson scheme given by X n+1 -X n-1 2∆t = 1 6 (A n+1 + 4 A n + A n-1 ).


Jose Ordonez-Miranda

Laurent Tranchant

Sergei Gluchko

Thomas Antoni

Sebastian Volz
email: sebastian.volz@ecp.fr

Fresnel-like Formulas for the Reflection and Transmission of Surface Phonon-Polaritons at a Dielectric Interface

The reflection and transmission coefficients of a surface phonon-polariton propagating along the surface of a thin film of SiO 2 and crossing the interface of two dielectric media are analytically determined. Based on the expansion of the electrical and magnetic fields in terms of normal modes, explicit expressions for the reflectivity and transmissivity of the radiation fields generated at the dielectric interface are also obtained. Symmetrical and simple Fresnel-like formulas are derived for nanofilms. For the dielectric interfaces of air/BaF 2 and air/Al 2 O 3 , it is shown that: i) The polariton reflectivity (transmissivity) decreases (increases) as the film thickness increases, while its radiation equivalent follows the opposite behavior. ii)

In the polariton and radiation fields, the transmissivity is significantly more sensitive than the reflectivity to the changes on the permittivity mismatch of the dielectric interface. For a 143 nm-thick film, the polariton transmissivity (reflectivity) changes 13.2% (1.9%), when this mismatch varies by 50%. iii) The reflectivity and transmissivity of the radiation fields are smaller than their polariton counterparts, which together account for around 82% of the total energy. The proposed formalism accurately fulfils the principle of conservation of energy for describing the reflection and transmission of both the polariton and radiation fields generated at a dielectric interface.

I. INTRODUCTION

Surface phonon-polartions (SPPs) have been recently proposed [1][2][3] as potential energy carriers to enhance significantly the phonon heat transport in polar nanomaterials with low thermal conductivity, as is the case of the amorphous SiO 2 and crystalline SiC. Many theoretical [4][5][6][7] and experimental [8,9] studies have shown that the propagation length of SPPs traveling along the interface of nanofilms can be more than three orders of magnitude larger than the corresponding phonon mean free path, which yields in turn a SPP thermal conductivity comparable or higher than the phonon counterpart. We have also shown that the thermal conductance of polar nanowires due to SPPs is quantized and its value can be comparable to the one of phonons at room temperature [10]. Furthermore, various research groups have shown that the SPPs have promising applications on the thermal performance of nanoscale devices [4,8,11], radiative heat transfer [12][13][14][15], high-density infrared data storage [16], surface infrared absorption [17], coherent thermal emission [18], and photonics [19,20]. Given that SPPs are electromagnetic waves that propagates along the interface between polar and dielectric materials, their energy transport is expected to increase as the material size is scaled down to nanoscales, due to the predominance of the surface effects at these size scales [4,11]. SPP energy transport is determined by material permittivities and hence it can be modified by material discontinuities, which can be present in optical/electronic integrated circuits printed out on surfaces supporting the propagation of SPPs. The study of the reflection and transmission of SPPs at these discontinuities is thus of great practical interest and it has attracted considerable experimental [21][22][23] and theoretical [24][25][26] interest in recent years. By using Fourier transform infrared spectroscopy, Chen et al. [21] measured the transmissivity through a 1 µm-thick film of SiO 2 deposited over a substrate of a 500 µm-thick silicon wafer and showed that the minima of the fraction of transmitted energy occur at the resonance wavelengths of absorption. For the major resonance peak, they found a fivefold increase in the transmission across a solid film through a perforated film with a periodic array of circular holes. These results are consistent with the numerical ones obtained by Catrysse and Fan [24], for the transmission through subwavelength hole arrays in a thin film of SiC, by means of a three-dimensional 3D finite-difference time domain method. Similar experimental data were reported by Gall et al. [22] and Marquier et al. [23] for the reflectivity and transmissivity spectra in a SiC grating supporting the propagation of SPPs. Furthermore, by placing infinitely conducting metallic planes both above and below the propagation interface and making a waveguide-like structure, Stegeman et al. [25][26][27] numerically calculated the transmission and reflection coefficients for a surface plasmon polariton crossing a plane boundary at normal incidence. They showed that these coefficients vary strongly with the ratio between the separation distance of the metallic planes and the polariton wavelength, such that the reflection (transmission) coefficient tends to zero (unity) as this ratio increases.

Even though the formalism proposed by these latter authors is rigorous and could also be extended for SPPs, their numerical calculations are quite complicated and time consuming to implement, as was mentioned by the same authors [27].

In this paper, we analytically determine the reflection and transmission coefficients of a SPP propagating along the surface of a thin film of SiO 2 and crossing the interface of two dielectric media. Based on the expansion of the electrical and magnetic fields in terms of normal modes, explicit expressions for the reflection and transmission coefficients of the radiation fields generated at the dielectric interface are also derived. It is shown that these latter fields are indispensable for satisfying the boundary conditions and the principle of conservation of energy. Our results provide clear physical insights on how the SPP and radiation energy depends on the permittivity of the dielectric media as well as on the permittivity and thickness of the thin film.

II. FORMULATION OF THE PROBLEM AND NORMAL MODES

Let us consider a layered system supporting the propagation of a SPP, as shown in Fig. 1.

The thin film of thickness d and permittivity ε 1 (ω) is a polar material (as SiO 2 and SiC) deposited over a metallic substrate, which limits the generation of radiation modes at the vertical interface of the semi-infinite dielectric media of permittivities ε 2 and ε 3 [28]. When the incident SPP propagating along the interface 1-2 arrives to the vertical interface z = 0, it is partially reflected backward and transmitted to the interface 1-3. In general, the incident SPP can also be partially reflected and transmitted as electromagnetic radiation, which is predicted by the Maxwell equations and is required to satisfy the boundary conditions at z = 0, as shown below. The problem consists in finding the reflection and transmission coefficients of both the SPP and radiation electromagnetic fields, which will allow us to quantify the distribution of energy in the reflected and transmitted SPP, as well as in the radiation losses. By analogy with the problems involving waveguides, the electrical and magnetic fields for the SPP and radiation fields in each region of the system shown in Fig. 1 can be expressed in terms of the corresponding modes in absence of the discontinuity at z = 0 (ε 2 = ε 3 ). By solving the Maxwell equations under the transverse magnetic polarization, which is required for the existence of the incident SPP [4], it can be shown that the non-zero components of these fields are given by [START_REF] Yeh | The Essence of Dielectric Waveguides[END_REF]

E (n) x = β ωε H (n) y , (1a)
E (n) z = 1 iωε n ∂H (n) y ∂x , (1b)
∂ 2 H (n) y ∂x 2 -p 2 n H (n) y = 0, ( 1c
)
where β and p n are the in-plane and transverse wave vectors, respectively, and they are related by

p 2 n = β 2 -ε n k 2 0 , k 0 = ω/c
, ω is the excitation frequency, c is the speed of light in vacuum, ε n is the permittivity of the medium n = 1, 2, and i = √ -1 is the imaginary unit. The electrical

E (n) x,z
and magnetic H

(n) y components depend on the coordinate z and time t through their common factor exp(i(ωt -βz)), which is omitted in Eqs. (1a)-(1c)

and will be suppressed hereafter. For ε 3 = ε 2 , Eq. (1c) establishes that the magnetic field inside the film H

(1) y and the semi-infinite medium H

(2) y are given by

H (1) y = A cosh (p 1 (x + d)) , (2a)
H (2) y = Be p 2 x + Ce -p 2 x , (2b)
where H

(1) y has been chosen in such a way that the associated electrical field E

x satisfies the boundary condition E x (x = -d) = 0, which is imposed by the metallic substrate [START_REF] Jackson | Classical electrodynamics[END_REF]. The parameters A, B, and C are determined by the excitation source and the boundary conditions H

(1)

y = H (2) y and E (1) z = E (2) z
at the interface x = 0. These latter relations yield

B + C = A cosh(p 1 d), (3a)
B -C = Aα 21 sinh(p 1 d), (3b)
where α 21 = ε 2 p 1 /ε 1 p 2 . Equations ( 1)-( 3) are valid for both the SPP and radiation fields, which are going to be analyzed separately, as follows.

A. SPP modes

These modes are defined by the exponential decay (Re(p 2 ) > 0) of the fields as they travel away from the interface x = 0. This condition requires that B = 0, for which Eqs. (3a) and (3b) yield the following dispersion relation

p 2 ε 2 + p 1 ε 1 tanh(p 1 d) = 0. (4)
According to Eqs. ( 2) and (3a), the transverse magnetic field (h(x) = H y ) can be conveniently written as

h(x) = N S    cosh(p 1 (x + d)), -d < x < 0 cosh(p 1 d)e -p 2 x , x > 0 (5)
where N S = A for the SPP modes, and the transverse electrical field is given by e(x) = (β/ωε)h(x) (see Eq.(1a)). Equation ( 5) thus defines one SPP mode for each pair (p 1 , p 2 ) of transverse wave vectors that satisfy the dispersion relation in Eq. ( 4). According to the appendix A, these modes are orthogonal and they satisfy the normalization condition

∞ -d e(x)h * (x)dx = 1, (6a)
N 2 S = 2ωε 1 βdγ , (6b)
γ = 1 + sinh(2p 1 d) 2p 1 d + ε 1 ε 2 cosh 2 (p 1 d) p 2 d , (6c)
where the * stands for the complex conjugate of h(x).

B. Radiation modes

The electrical and magnetic fields of these modes are not spatially damped along the x direction and therefore they are defined by Re(p 2 ) = 0. In this case, the transverse wave vector p 2 can be defined as p 2 = -ik, where k > 0 is the radiation wave vector along the transverse x direction. Based on Eqs. ( 2) and (3) and by analogy with the SPP modes, the radiation modes can be defined as follows

h(x, k) = N R    cosh(p 1 (x + d)) , -d < x < 0 cosh(p 1 d) cos(kx) + ξ sinh(p 1 d) sin(kx) , x > 0 (7a)
where N R = A for the radiation modes, ξ = ε 2 p 1 /ε 1 k, and

β 2 = ε 2 k 2 0 -k 2 = p 2 1 + ε 1 k 2 0 .
The corresponding electrical field is given by e(x, k) = (β/ωε)h(x, k), as before. These radiation modes are also orthogonal and they are normalized by the following condition (Appendix

A) ∞ -d e(x, k)h * (x, k )dx = δ(k -k ), (8a)
N 2 R = 2ωε 2 πβ(k)γ R , ( 8b
) γ R = | cosh(p 1 (k)d)| 2 + |ξ sinh(p 1 (k)d)| 2 . ( 8c
)
We can now use the previous SPP and radiation modes to express the electrical and magnetic fields involved in the system shown in Fig. 1. Let e < 0 (x) and h < 0 (x) (e > 0 (x) and h > 0 (x)) be the transverse components of the electrical and magnetic fields of the incident (transmitted) SPP, respectively; and e < (x, k) and h < (x, k) (e > (x, k) and h > (x, k)) be the corresponding electromagnetic fields of the backward (forward) radiation modes. The transverse components of the total electromagnetic fields in the regions z < 0 (E < x , H < y ) and z > 0 (E >

x , H > y ) can then be expressed as follows [START_REF] Shevchenko | Continuous transitions in open waveguides[END_REF]]

E < x (x, z) = e < 0 (x) e -izβ < -re izβ < - ∞ 0 r R (k)e < (x, k)e izβ < R dk, (9a)
H < y (x, z) = h < 0 (x) e -izβ < + re izβ < + ∞ 0 r R (k)h < (x, k)e izβ < R dk, (9b)
E > x (x, z) = τ e > 0 (x)e -izβ > + ∞ 0 τ R (k)e > (x, k)e -izβ > R dk, (9c)
H > y (x, z) = τ h > 0 (x)e -izβ < + ∞ 0 τ R (k)h > (x, k)e -izβ < R dk, (9d)
where r and τ (r R and τ R ) are the reflection and transmission coefficients of the SPP (radiation) fields, respectively, and β < and β > (β < R and β > R ) are the SPP (radiation) wave vectors along the propagation direction, for z < 0 and z > 0, respectively. The continuity of the

tangential electrical (E < x (x, 0) = E > x (x, 0)) and magnetic (H < y (x, 0) = H > y (x, 0)) fields at the dielectric interface z = 0 yields (1 -r) e < 0 (x) - ∞ 0 r R (k)e < (x, k)dk = τ e > 0 (x) + ∞ 0 τ R (k)e > (x, k)dk, (10a)
(1 + r) h < 0 (x) + ∞ 0 r R (k)h < (x, k)dk = τ h > 0 (x) + ∞ 0 τ R (k)h > (x, k)dk. (10b)
Note that in absence of the radiation terms, Eqs. (10a) and (10b) cannot be satisfied because of the transverse wave vectors in the region z < 0 are different than the corresponding ones in z > 0. It is therefore clear that the presence of the radiation modes are indispensable to satisfy these boundary conditions. This mode matching approach was previously applied to numerically determine the reflection and transmission coefficients of surface plasmon-polaritons at dielectric interfaces [START_REF] Oulton | [END_REF]33]. To eliminate the dependence on the variable x and to solve Eqs. (10a) and (10b) for the coefficients r, τ , r R , and τ R , it is convenient to define the following inner product

(f, g) = ∞ -d f (x)g * (x)dx, (11)
where f and g are functions of x over the interval of integration. Multiplying Eq. ( 10a) by (h < 0 ) * and Eq. ( 10b) by (e > 0 ) * , and integrating both sides over the interval -d < x < ∞;

the following expressions are found for τ

τ = (1 -r) (e < 0 , h > 0 ) - ∞ 0 r R (k)(e < (k), h > 0 )dk, (12a)
τ = (1 + r) (h < 0 , e > 0 ) + ∞ 0 r R (k)(h < (k), e > 0 )dk, (12b)
where we have used the normalization condition in Eq.(6a) for the SPP modes and the orthogonality of the SPP and radiation modes, which is discussed in appendix A. By following a similar procedure, it is easy to show that τ R is determined by

τ R (k) = (1 -r) (e < 0 , h > (k)) - ∞ 0 r R (k )(e < (k ), h > (k))dk , (13a) τ R (k) = (1 + r) (h < 0 , e > (k)) + ∞ 0 r R (k )(h < (k ), e > (k))dk . (13b)
Equations ( 12) and (13) show that the exact determination of the reflection and transmission coefficients is quite complicated due to the presence of the integral terms. However, for dielectric media with permittivities ε 2 and ε 3 of the same order of magnitude, both reflection coefficients are expected to be close to zero (r, r R → 0), and the electromagnetic modes in the regions z < 0 and z > 0 are approximately orthogonal

[(e < (k), h > 0 ), (h < (k), e > 0 ), (e < 0 , h > (k)), (h < 0 , e > (k))] → 0 and (e < (k ), h > (k)) ≈ δ(k -k ).
Under these conditions and a first-order approximation of Eqs. ( 12) and ( 13), it is found that the reflection and transmission coefficients of the SPP and radiation modes are analytically given by

r = (e < 0 , h > 0 ) -(h < 0 , e > 0 ) (e < 0 , h > 0 ) + (h < 0 , e > 0 ) , (14a)
τ = 2(e < 0 , h > 0 )(h < 0 , e > 0 ) (e < 0 , h > 0 ) + (h < 0 , e > 0 ) , ( 14b
) 2r R (k) = (e < 0 , h > (k)) -(h < 0 , e > (k)), (14c) 2τ R (k) = (e < 0 , h > (k)) + (h < 0 , e > (k)). ( 14d
)
Note that while the SPP reflection and transmission coefficients depend only on inner products of SPP modes, their radiation counterparts are determined by inner products between the SPP modes in z < 0 with radiation modes in z > 0. This means that the origin of the radiation fields is the incident SPP field, as expected. After evaluating the inner products in Eqs. (14a)-(14d) with the help of Eqs. ( 5)-( 8) and ( 11), the reflection and transmission coefficients can be written explicitly as follows

r = ε 2 ε 3 (β < -β > )I 1 + ε 1 (ε 3 β < -ε 2 β > )I 2 ε 2 ε 3 (β < + β > )I 1 + ε 1 (ε 3 β < + ε 2 β > )I 2 , ( 15a
) τ = σ (ε 2 I 1 + ε 1 I 2 )(ε 3 I 1 + ε 1 I 2 ) ε 2 ε 3 (β < + β > )I 1 + ε 1 (ε 3 β < + ε 2 β > )I 2 , (15b) r R = σ (β < -(β > R ) * ) J 1 ε 1 + β < ε 2 - (β > R ) * ε 3 J 2 , ( 15c
) τ R = σ (β < + (β > R ) * ) J 1 ε 1 + β < ε 2 + (β > R ) * ε 3 J 2 , (15d)
where the parameters σ = 4

β < β > /δ < δ > , σ R = ε 1 ε 3 d/(β < (β > R ) * δ < η > )
, and the dimentionless integrals

I 1 = 1 2 sinh [(p < 1 + p > 1 )d] (p < 1 + p > 1 )d + sinh [(p < 1 -p > 1 )d] (p < 1 -p > 1 )d , (16a)
I 2 = cosh(p < 1 d) cosh(p > 1 d) (p 2 + p 3 )d , (16b)
J 1 = 1 2 sinh [(p < 1 + (p > 1R ) * )d] (p < 1 + (p > 1R ) * )d + sinh [(p < 1 -(p > 1R ) * )d] (p < 1 -(p > 1R ) * )d , (16c)
J 2 = cosh(p < 1 d) (p 2 2 + k 2 )d p 2 cosh((p > 1R ) * d) + ε 3 ε 1 (p > 1R ) * sinh((p > 1R ) * d) . (16d)
Equations ( 15a)-(15d) clearly show that the reflection and transmission coefficients for the SPP and radiation fields depend strongly on the permittivity mismatch of the semi-infinite media, as well as on the permittivity and thickness of the thin film. In absence of the dielectric interface (ε 3 = ε 2 ), the in-plane wave vectors β < = β > and therefore r = 0, τ = 1, and r R = τ R = 0, as expected. For very thin (k 0 d << 1) and very thick (k 0 d >> 1) films, Eqs. (15a) and (15b) becomes independent of the film thickness d and they can be significantly simplified. In the first limiting case, Eqs. (15a) and (15b) reduce to

r 0 = √ ε 3 - √ ε 2 √ ε 3 + √ ε 2 , ( 17a
)
τ 0 = 4 4 √ ε 2 ε 3 √ ε 2 + √ ε 3 ε 2 ε 3 (ε 2 -ε 1 )(ε 3 -ε 1 ) ε 2 (ε 2 -ε 1 ) + ε 3 (ε 3 -ε 1 ) , ( 17b
)
where the subscript "0" indicates that the formulas are valid for (k 0 d << 1). For thick films, on the other hand, Eqs.(15a) and (15b) take the form

r ∞ = ε 2 ε 3 (β < -β > ) + ε 1 (ε 3 β < -ε 2 β > )χ ε 2 ε 3 (β < + β > ) + ε 1 (ε 3 β < + ε 2 β > )χ , ( 18a
)
τ ∞ = ψ (ε 2 + ε 1 χ)(ε 3 + ε 1 χ) ε 2 ε 3 (β < + β > ) + ε 1 (ε 3 β < + ε 2 β > )χ , (18b)
χ = - ε 1 (p < 1 + p > 1 ) ε 2 p < 1 + ε 3 p > 1 , (18c)
ψ = - 4ε 2 ε 3 β < β > p < 1 p > 1 (p < 1 + p > 1 ) (ε 2 1 -ε 2 2 )(ε 2 1 -ε 2 3 ) , ( 18d
)
where

β < = k 0 ε 1 ε 2 /(ε 1 + ε 2 ), β > = k 0 ε 1 ε 3 /(ε 1 + ε 3 ), p < 1 = -ε 1 /ε 2 β < , p > 1 = -ε 1 /ε 3 β > , and the subscript "∞" stands for k 0 d >> 1. Note that both r 0 (ε 2 , ε 3 ) = -r 0 (ε 3 , ε 2 ) and r ∞ (ε 2 , ε 3 ) = -r ∞ (ε 3 , ε 2 )
, which means that when the semi-infinite media are interchanged, the reflection coefficient of the reflected SPP changes by a phase reversal only. This behavior coincides with the one exhibited by the corresponding Fresnel equation for normal incidence of plane waves [START_REF] Jackson | Classical electrodynamics[END_REF]. By contrast, this phase reversal is not present in the transmission coefficient, which is invariant under the interchange of the semi-infinite

media (τ 0 (ε 2 , ε 3 ) = τ 0 (ε 3 , ε 2 ) and τ ∞ (ε 2 , ε 3 ) = τ ∞ (ε 3 , ε 2 )).
Given that the energy fractions of the reflected and transmitted SPP are given by the square of the magnitude of r and τ , respectively; as shown below, the aforementioned features of the SPP reflection and transmission coefficients establish that the energy of the reflected and transmitted SPP keeps invariant under the interchange of the dielectric media. This feature is exclusive of SPPs and is not present in the radiation fields, as indicated by Eqs. (15c) and (15d).

The accuracy of the simple and approximate expressions for the reflection and transmission coefficients in Eqs. (14a)-(14c) can be evaluated by calculating the integrals in Eqs. ( 12) and ( 13) with the estimated value of the refleccion coefficient of the radiation modes in Eq. (14c). In doing this, the following expressions for the differences ∆r, ∆τ , and ∆τ R (k) between the exact expressions of the reflection and transmission coefficients (Eqs. ( 12) and ( 13)), and their corresponding approximate results (Eqs. (14a), (14b), and (14c));

are obtained ∆r = - ∞ 0 r R (k) [(e < (k), h > 0 ) + (h < (k), e > 0 )] dk (e < 0 , h > 0 ) + (h < 0 , e > 0 ) , (19a) ∆τ
= ∞ 0 r R (k) [(e < 0 , h > 0 )(h < (k), e > 0 ) -(h < 0 , e > 0 )(e < (k), h > 0 )] dk (e < 0 , h > 0 ) + (h < 0 , e > 0 ) , ( 19b
) ∆τ R (k) = 1 2 ∞ 0 r R (k ) [(h < (k ), e > (k)) -(e < (k ), h > (k))] dk . (19c)
Given the complicate dependence of r R (k) on the wave vector k (Eq. ( 15c)), the integrals in Eqs. (19a)-(19c) cannot be evaluated analytically. For very thin film (k 0 d << 1), which is the case of major interest to enhance the propagation of SPPs [4,11], these integrals simplify significantly and the deviations exhibit the following behavior ∆r and ∆τ

∝ (ε 2 -ε 1 )(ε 3 -ε 1 )(k 0 d) 3 (20a) ∆τ R (k) ∝ (ε 2 -ε 1 )(k 0 d) 2 , (20b)
which explicitly indicate that for k 0 d << 1, our approximate results in Eqs. (14a)-(14d) are accurate and their accuracy is expected to increase as the film thickness decreases. The fact that the SPP coefficients are proportional to the third power of the film thickness, while radiation one is proportional to the second power of this thickness; shows that the approximate expression of the reflection and transmission coefficients of the SPP modes become more accurate than the transmission coefficient of the radiation modes, as the film thickness is scaled down.

The power flow per unit film width of the reflected and transmitted SPP and radiation fields can now be calculated through the Poynting vector S along the propagation direction (z axis), which is defined by

S = 1 2 Re ∞ -d E x (x, 0)H * y (x, 0)dx , (21)
Based on Eqs. ( 9) and ( 21), we can determine the power of the incident (S inc ), reflected (S ref ), and transmitted (S tran ) SPP; as well as the power of the reflected (S ref R ) and transmitted (S tran R ) radiation fields. The reflectivity and transmissivity of the SPP and radiation fields are therefore defined as R = -S ref /S inc , T = S tran /S inc , R R = -S ref R /S inc , and T R = S tran R /S inc . In terms of the reflection and transmission coefficients, these relations yield

R = |r| 2 , (22a)
T = |τ | 2 , (22b) R R = ∞ 0 |r R (k)| 2 dk, (22c)
T R = ∞ 0 |τ R (k)| 2 dk. (22d)
Given that all the energy of the incident SPP is transfered to the reflected and transmitted SPP and radiation fields, the principle of conservation of energy establishes that

R + T + R R + T R = 1. ( 23
)
We will show below that the approximate, simple, and analytical results for the reflection and transmission coefficients involved in Eqs. (22a)-( 22d) are able to satisfy Eq. ( 23) with an error less than 2%.

III. RESULTS AND DISCUSSIONS

The reflectivity and transmissivity of the SPP and radiation fields generated by an incident SPP crossing the dielectric interface shown in Fig. 1 are quantified and analyzed in this section. Special emphasis is put on the effects of the film thickness and permittivities of the dielectric media. We consider a film of amorphous SiO 2 (silica), which is an abundant polar material in nature and widely used in electronic applications [1,4]. Furthermore, SiO 2 is commonly employed in the manufacturing of electrical circuits, in which thermal sources generated by Joule effect might be cooled via SPPs. For the dielectric media, we are going to use air, as a natural surrounding medium, in contact with BaF 2 or Al 2 O 3 (alumina), which are transparent crystals in a wide range of frequencies [START_REF] Palik | Handbook of Optical Constants of Solids[END_REF] and are used in many applications related to optics, infrared spectroscopy, and electronics. All of these materials behave as lossless materials (with real permittivities) at the frequency ω = 210 • 10 12 rad/s (k 0 = 0.7 µm -1 ), which is going to be used in this work, to guaranty the orthogonality of the SPP and radiation modes and to simplify the energy distribution analysis based on the Poynting vector, as discussed in appendix A. as shown in Fig. 1. Calculations were done using the dispersion relation in Eq. ( 4) and its analysis in Appendix B. It is seen that β inscreases with the film thickness, such that it reaches asymptotic values for very thin (k 0 d << 1) and very thick (k 0 d >> 1) films.

The thick-film limit is given by the well known dispersion relation

β = ε 1 ε 2 /(ε 1 + ε 2 )k 0 ,
and the thin-film one is described by the simple solution of the dispersion relation in Eq. [START_REF] Yeh | The Essence of Dielectric Waveguides[END_REF]. Given that k 0 = 0.7 µm -1 , this latter limit holds for films with thicknesses (d <<1.4 µm) from hundreds of nanometers downwards, which represent the regime of major energy transport by SPPs [6,7,11]. Note that irrespective of the the film thickness, the in-plane wave vector increases with the permittivity of the dielectric media and it is greater than the one in vacuum (β > k 0 ), which indicates that the SPP is propagating with a wave vector to the right of the light line [11]. Furthermore, it is clear that β → k 0 , when the thickness and/or the dielectric permittivity reduces, which implies a photon-like nature of the SPP.

As these parameters increase, β/k 0 separates from the unity and tends to a phonon-like behavior.

The film-thickness dependence of the normalized transverse wave vectors in the regions z < 0 (p < 1 , p 2 ) and z > 0 (p > 1 , p 3 ) are shown in Fig. 3, for the air/BaF 2 dielectric interface. and crossing a dielectric interface of air/BaF 2 (see Fig. 1), as a function of the normalized film thickness. Calculations were performed with the data used in Fig. 2.

As in the case of the in-plane wave vector β, the transverse wave vectors increase with the film thickness and they become independent of this thickness for very thin (k 0 d << 1) and very thick (k 0 d >> 1) films. For a given film thickness, these wave vectors are larger in the medium with higher absolute value of permittivity, as established by the dispersion relation in Eq. ( 4) and its analytical solution in Eq. ( 29). The fact that all transverse wave vectors are positive indicates that the electrical and magnetic fields of the SPPs decay spatially as they travel away from the propagation interface x = 0, which guaranties the existence of SPPs [4,11]. Given that |ε(SiO 2 )| > ε(BaF 2 )> ε(air), the transverse spatial attenuation inside the thin film is stronger than the one in BaF 2 , which in turns is more intense than the one within air. It is therefore clear that the SPP fields within the SiO 2 film are more confined to the propagation interface than those in the dielectric media. for micro-films, they are determined by a little more complicated but still analytical formulas in Eqs. (18) and (20). Equations ( 17) and ( 18) can thus be considered as the Fresnellike equations for the reflection and transmission of SPPs across the dielectric interface shown in Fig. 1. Note that the reflectivity (transmissivity) for the dielectric interface of Air/BaF 2 is higher (lower) than the corresponding one of the interface Air/Al 2 O 3 . This is quite reasonable due to the fact that the permittivity mismatch between the dielectric media in the Air/BaF 2 interface (ε 3 -ε 2 = 1) is higher than that for the another interface (ε 3 -ε 2 = 0.5). For this modification on the dielectric interface (∆ε 3 = 0.5) and k 0 d = 0.1 (d = 143 nm), the change on the reflectivity is ∆R = 1.9%, which rises to ∆T = 13.2%

for the transmissivity. The SPP transmissivity is therefore more sensitive than the SPP reflectivity to the changes on the permittivity mismatch of the dielectric interface. respectively, as a function of the normalized radiation wave vector k/k 0 . Both spectra follow a similar behavior, which does not vary significantly for film thicknesses k 0 d > 1 and k/k 0 < √ ε 3 . Within this range of the radiation wave vector, the wave vector β > R is real, as established just below Eq. (7a); while for k/k 0 > √ ε 3 , β > R becomes an imaginary number, which leads to the oscillations of the spectra. In this latter interval of the radiation wave vector, both spectra decay fast to zero as k/k 0 increases; which indicates that the major contributions to R R and T R (area under the curves) arise from wave vectors within the interval k/k 0 < √ ε 3 .

1 0 -2 1 0 -1 1 0 0 1 0 1 1 0 -5 1 0 -4 1 0 -3 1 0 -2 1 0 -1 A i r / B a F 2
k 0 d = 1 k 0 d = 5 √ 3 1 0 -2 1 0 -1 1 0 0 1 0 1 1 0 -4 1 0 -3 1 0 -2 1 0 -1 A i r / B a F 2
Note that, for a given film thickness, the values of the transmissivity spectrum are higher than those of the reflectivity one, and hence the transmissivity is espected to be higher than the reflectivity (T R > R R ). This is confirmed by Fig. 7, not only for the Air/BaF 2 dielectric interface but also for the Air/Al 2 O 3 one. For both interfaces, the radiation reflectivity Figure 8 shows the film-thickness dependence of the reflected and transmitted power fractions involved in the SPP (R + T ) and radiation (R R + T R ) fields, as well as their overall sum, for the Air/BaF 2 dielectric interface. The fact that R + T >> R R + T R means that the major part of energy of the incident SPP is distributed between the reflected and trasmitted SPPs, while less than 20% of its energy goes to the radiation fields. Note that R + T + R R + T R ≈ 1, which satisfies quite well the principle of conservation of energy, with a deviation of less than 2% that occurs at k 0 d = 0.36 (d=514 nm). This deviation reduces to less than 0.5% for the dielectric interface of Air/Al 2 O 3 , due to its smaller permittivity mismatch than the one of the Air/BaF 2 interface. It is therefore clear that, when the permittivities of the dielectric media are of the same order of magnitude, the proposed formalism is able to analytically describe with good accuracy the reflection and transmission of both the SPP and radiation fields generated at a dielectric interface. Given that Eqs.

(15) and ( 20) have been derived assuming that the permittivivities of the dielectric media are not so dissimilar, in general, this simple approach is expected to yield deviations higher than 2%, from the conservation of energy, as the difference of permittivities of the dielectric media increases to values larger than the unity (ε 3 -ε 2 > 1). In these cases, Eqs. (15a)-(15d)

are not longer suitable and the determination of the reflection and transmission coefficients is much more complicated due to the presence of the radiation modes, as established by Eqs. ( 12) and ( 13), which still hold and can be conveniently solved by means of a numerical approach.

IV. CONCLUSIONS

The reflection and transmission of a surface phonon-polariton propagating along the surface of a thin film of SiO 2 and crossing the interface of two dielectric media has been analysed by means of an analytical approach based on the expansion of the electrical and magnetic fields in terms of normal modes. Fresnel-like formulas for the reflectivity and transmissivity of both the polariton and radiation fields generated at the dielectric interface have been explicitly determined. For the dielectric interfaces of air/BaF 2 and air/Al 2 O 3 , it has been shown that: i) The polariton reflectivity (transmissivity) decreases (increases) as the film thickness increases, while its radiation equivalent follows the opposite behavior. ii)

In the polariton and radiation fields, the transmissivity is significantly more sensitive than the reflectivity to the changes on the permittivity mismatch of the dielectric interface. iii)

The reflectivity and transmissivity of the radiation fields are smaller than their polariton counterparts, which together account for around 82% of the energy of the incident surface phonon-polariton. The proposed formalism accurately fulfils the principle of conservation of energy and could be used for quantifying the polariton energy and radiation losses at a dielectric interface.

The right-hand side of Eq.( 24) can be evaluated and simplified by splitting it for positions inside the thin film (-d < x < 0) and the dielectric medium (x > 0), applying the corresponding boundary conditions at x = -d and x = 0, and taking into account the disappearance of the fields at x = ∞. The final result is

p 2 n -(p 2 m ) * iβ n ∞ -d e n h * m dx = ε * 1 ε 1 - ε * 2 ε 2 e (1) * zm h (1) n | x=-d . (27)
Equation ( 25) thus indicates that when the permittivities of the thin film and the dielectric medium are real, the integral in its left-hand side vanishes

) 28
for any n = m. This orthogonality relation holds for any combination of SPP modes as well as for an inner product between SPP and radiation modes. Taking into account the definition of the Poynting vector in Eq.( 19), Eq.( 26) establishes that the power cannot be

FIG. 1 .

 1 FIG.1. SPP and radiation modes generated by an incident SPP at a dielectric interface. The metallic substrate of the film of thickness d is used to diminish the radiation modes and enhance the SPP energy transport.

Figure 2 3 FIG. 2 .

 232 Figure2shows the film-thickness dependence of the wave vector β along the propagation direction of a SPP traveling via the interface of a film of SiO 2 and a dielectric medium,

3 T r a n s v e r s e w a v e v e c t o r s , p n / k 0 N

 30 FIG. 3. Normalized transverse wave vector of a SPP propagating along the surface of a SiO 2 film

Figure 4 FIG. 4 .

 44 Figure4shows that reflectivity R and transmissivity T of the SPP fields as a function of the normalized film thickness and for two dielectric interfaces. For both cases, R decreases as the film thickness increases, while T exhibits the oposite trend; and they tend to thickness independent values for very thin (k 0 d << 1) and very thick (k 0 d >> 1) films, which is consistent with the behavior of the SPP wave vectors shown in Figs.2 and 3. For nanofilms, R and T are described by the simple and symmetrical formulas in Eqs. (17) and (20), while

  FIG. 5. Deviations of the reflection and transmission coefficients of the SPP modes propagating along the surface of a SiO 2 film and crossing a dielectric interface, as a function of the normalized film thickness. Calculations were done for two dielectric interfaces of Air/BaF 2 and Air/Al 2 O 3 , and with the data used in Fig. 2.

R a d i a t i o n r e f l e c t i v i t y s p e c t r a , |r ( k ) | 2 N

 2 o r m a l i z e d w a v e v e c t o r , k / k 0 k 0 d = 1/5

R a d i a t i o n t r a n s m i s s i v i t y s p e c t r a , |ؔ ( k ) | 2 N o r m a l i z e d w a v e v e c t o r , k / k 0 k 0 d = 1/ 5 k 0 d = 1 k 0 d = 5 √ 3 FIG. 6 .

 2051536 FIG. 6. Spectra of the (a) radiation reflectivity and (b) radiation transmissivity as a function of the normalized radiation wave vector. Calculations were performed for three thicknesses of a film of SiO 2 and a dielectric interface of Air/BaF 2 .

(

  transmissivity) increases (decreases) as the film thickness increases, which is opposite to the behavior exhibited by the SPP counterpart shown in Fig.4. This is reasonable given that the energy of the incident SPP must be conserved. Both R R and T R increase as the permittivity mismatch of the dielectric interface increases from ε 3 -ε 2 = 1/2, (for Air/Al 2 O 3 ) to ε 3 -ε 2 = 1 (for Air/BaF 2 ), as expected. For this change of the dielectric interface (∆ε 3 = 0.5) and k 0 d = 0.1 (d = 143 nm), the variation of the reflectivity is ∆R R = 0.5%, which rises to ∆T R = 12.2% for the transmissivity. These changes are smaller than the corresponding ones of the SPP reflectivity and SPP transmissivity shown in Fig.4. Furthermore, the comparison of Figs.4 and 6, for a particular thickness and dielectric interface, shows that R > R R and T > T R , which indicates that reflected and transmitted SPP fields carry more energy than their radiation equivalents.

2 AFIG. 7 . 2 FIG. 8 .

 2728 FIG. 7. Reflectivity and Transmissivity of the radiation fields generated by a SPP propagating along the surface of a SiO 2 film and crossing a dielectric interface, as a function of the normalized film thickness. Calculations were done for two dielectric interfaces of Air/BaF 2 and Air/Al 2 O 3 , and with the data used in Fig. 2.

APPENDIXA.

  Orthogonality of modesLet us consider that (e n , h n ) and (e m , h m ) are two different modes propagating along the system shown in Fig.1with ε 3 = ε 2 , they then satisfy Eqs. (1a) and (1c)∂ 2 h n ∂x 2 -p 2 n h n = 0,(24a)∂ 2 h m ∂x 2 -p 2 m h m = 0. (24b)We can see that Eqs.(22a) and (22b) define a classical Sturm-Liouville problem, and hence the orthogonality of the n th and m th modes can be easily proved multiplying Eq. (22a) and the complex conjugate of Eq. (22b) by h * m and h n , respectively. The subtraction of the resulting equations yieldsp 2 n -(p 2 m ) * h n h * m = . (1a) (e = (β/ωε)h) and (1b) (e z = (1/iωε)∂h/∂x) for each mode and integrating both sides of Eq.(23), we obtain the following

transferred among different modes in lossless media. Even though real materials always exhibit losses, these becomes negligible at certain excitation frequencies, as is the case of SiO2 and SiC [11,[START_REF] Palik | Handbook of Optical Constants of Solids[END_REF]. This is the reason why in this work, we are considering lossless media, for which Eq. ( 26) is valid. On the other hand, the normalization conditions of the SPP (Eq. (6a)) and radiation (Eq. (8a)) modes are straightforward to determine by direct integration of the modes in Eqs. ( 5) and (7), respectively.

B. Solution of the SPP dispersion relation

By writing the in-plane wave vector as β = √ εk 0 , the transverse wave vectors p n = √ ε -ε n k 0 and the SPP dispersion relation in Eq. ( 4) reads

where λ = k 0 d is the normalized film thickness. Considering that λ < 1, which is the case of interest to enhance the SPP energy transport [11], Eq.( 27) can be solved through perturbation theory. For an approximation up to λ 4 , Eq. ( 27) takes the form

which indicates that the effective permittivity ε of the system has the following expansion

The combination of Eqs. ( 29) and (30) yields the following values for the parameters ε (2) and ε (4)

Equations ( 24) and ( 25) are valid for any complex permittivity ε 1 of the film. For a lossless film (ε 1 is a real number), however, these equations hold for ε 1 < 0 only. This is the condition of existence and propagation of SPPs along very thin films (λ << 1), which differs from the well known constraint ε 1 < -ε 2 , valid for thick films (λ >> 1).


Jose Ordonez-Miranda

Laurent Tranchant

Beomjoon Kim

Yann Chalopin

Thomas Antoni

Sebastian Volz

Anisotropy and size effects on the thermal conductivity of polar nano thin films due to surface phonon-polaritons

The effects of the permittivity anisotropy and size of a nano thin film on its thermal conductivity due to surface phonon-polaritons are studied. We demonstrate that this thermal conductivity is a linear combination of the inverse first and third powers of the film thickness. For a 100 nm-thick film of α-quartz surrounded by air, it is shown that its thermal conductivity along the optical axis is comparable to the phonon counterpart and equals to 13 W/m.K, which is 25% higher than that along a perpendicular direction, at room temperature. Higher values are found for thinner films at higher temperatures.

Surface phonon-polaritons (SPPs) are electromagnetic waves generated by the coupling between photons and phonons at the interface between two media [1][2][3]. Over the past few years, various research groups have shown that these surface waves have promising applications on the thermal performance of nanoscale devices [4][5][6], radiative heat transfer [7][8][9][10], high-density infrared data storage [11], surface infrared absorption, [12] coherent thermal emission, [13] and photonics [14,15]. In these cases, the operating principles are based on the fact that in nanomaterials, the surface effects predominate over the volumetric ones, and therefore the energy transport by SPPs is particularly important. The SPP energy contribution increases as the material size is scaled-down in the direction perpendicular to the propagation one [4,6].

The nature of SPPs propagating along the interface of isotropic media is now well understood. SPPs can have long propagation lengths, of the order of a few centimeters on thin films [4][5][6], which in turn could lead to an efficient energy transport. Based on the Boltzmann transport equation, Chen et al. [4] showed that the SPP thermal conductivity of a 40 nm-thick thin film of silicon dioxide suspended in air, is 4 W/m.K at 500 K. This value represents an increase of about 100% over its bulk phonon thermal conductivity. This result has been confirmed by the fluctuation-dissipation theorem [16] and is consistent with experimental data for the propagation length [5]. More recently and under a fully analytical approach, we reconfirmed previous numerical results [4], for the SPP propagation length and thermal conductivity [6], and showed that a substrate with a permittivity slightly different than that of the superstrate of a nano thin film can generate large propagation lengths and therefore en- * Corresponding author: sebastian.volz@ecp.fr hance significantly its SPP thermal conductivity. Given that the absorption of energy by the film and its surrounding media is determined by their permittivities, it is expected that thin films with anisotropic permittivity can provide more variables than isotropic ones to control and optimize the SPP energy transport.

The purpose of this letter is to analyze the effects of the anisotropic permittivity and size of biaxial crystalline thin films on their SPP thermal conductivity. For nano thin films, an explicit expression for the SPP thermal conductivity as a function of their thickness is obtained, which provides clear physical insights on how the heat transport depends on the permittivity of the thin film and its surrounding media. The SPP propagation along the principal axes of the crystal are discussed in detail.

Let us consider a layered system supporting the propagation of SPPs, as shown in Fig. 1. The thin film is a biaxial crystal described by its permittivity tensor

ε =   ε x 0 0 0 ε y 0 0 0 ε z   , (1)
where ε j for j = x, y, z are the frequency-dependent permittivities along the principal axes of the anisotropic crystal. The surrounding medium on either side has a constant permittivity ε 1 > 0, as is the case of air, within a wide range of frequency [START_REF] Palik | Handbook of Optical Constants of Solids[END_REF]. By solving the Maxwell equations under proper bound-ary conditions for the transverse magnetic polarization required for the existence of surface polaritons [4], the following dispersion relation for the wave vector β, along the x axis is obtained

tanh (p 2 d/2) = - ε 1 p 2 ε x p 1 , (2)
where the transverse wave vectors inside and outside the film are given by p

2 2 = (ε x /ε z ) β 2 -ε z k 2 0 and p 2 1 = β 2 -ε 1 k 2 0 , respectively; k 0 = ω/c
, ω is the excitation frequency and c is the speed of light in vacuum. Note that the dispersion relation in Eq. ( 2) depends on the permittivity along (ε x ) and across (ε z ) the direction of SPP propagation, but not on the permittivity along the y axis, because the SPPs do not propagate in this direction. Equation (2) can significantly be simplified for a thin film (| p 2 | d << 2), which is of interest in this work to enhance the SPP propagation along the film. For a nano thin film of α-quartz, this condition is well satisfied for d ≤ 300 nm. In this case, Eq. (2) reduces to

p 2 2 ε x + 2 d p 1 ε 1 = 0, (3)
which, given the definition of p 2 , depends only on the permittivity ε z along the axis perpendicular to the film plane. As the film thickness is scaled down, the coupling of the electric and magnetic fields traveling along the two interfaces of the film becomes stronger, and therefore the dispersion relations of these fields should become more sensitive to ε z than to ε x . This shows that the separation of the permittivities in the dispersion relation of nano thin films is reasonable. Equation ( 3) is equal to our corresponding previous result for isotropic media after the substitution ε z → ε 2 [6]. Therefore, based on the analytical expressions for the real part β R of β = β R +iβ I and propagation length Λ = 1/ (2β I ) derived for this case [6], the SPP thermal conductivity can be written as

κ = A d 3 + B d , (4)
where

A = 1 2π ∞ 0 ω ξk 2 0 |ε z | 2 ε 2 1 ε I ∂f ∂T dω, (5a)
B = 1 8π ∞ 0 ω ξ ε 1 ε I - ε 1 ε I ξ ∂f ∂T dω, (5b)
being

ε z = ε R + iε I and ξ = |ε z | 2 -ε 1 ε R ,
which must be positive for any ω > 0, to guarantee the propagation of SPPs (Λ > 0). Within the thin film limit (| p 2 | d << 2), the simple Eq.( 4) establishes that the increase of the SPP thermal conductivity as the film thickness decreases is totally determined by the material-dependent parameters A and B. When 3 √ A << d < B, κ is driven by the inverse of the film thickness; while for B << d < 3 √ A, κ is inversely proportional to the cube of this thickness. The smaller values of the film thickness d are set by the surface roughness height and peak-to-peak distance, which are typically of the order of 1 nm. [4] Note that both parameters A and B, and hence the SPP thermal conductivity, increase as the permittivity ε 1 of the surrounding media decreases [4]. This result agrees with previous numerical simulations and is due to the fact that the propagation length Λ increases with ε 1 [6].

The quantification of the SPP thermal conductivity and its dependence on the anisotropic permittivity and size are now analyzed for a thin film of α-quartz surrounded by air (ε 1 = 1). This material is crystalline silica, which is common in nature and applications [1,[START_REF] Lamprecht | [END_REF]. The permittivities in the direction parallel (j = ) and perpendicular (j =⊥) to its optical axis are given by [START_REF] Lamprecht | [END_REF]

ε j (ω) = ε ∞ j Nj n=1 Ω 2 jn -ω 2 -iγ jn ω ω 2 jn -ω 2 -iγ jn ω , (6)
where Ω jn and ω jn are the longitudinal and transverse optical frequencies, respectively; γ jn is a damping constant and ε ∞ j is the high frequency permittivity. The experimental values of these parameters for the N = 4 (N ⊥ = 8) vibrational branches parallel (perpendicular) to the optical axis of α-quartz are reported in the literature [START_REF] Lamprecht | [END_REF]19]. Based on Eq. ( 6), it can be shown that the maxima of the imaginary part of ε j occurs at ω = ω jn , for ω jn >> γ jn , as is the case of α-quartz and most crystals. This indicates that these materials absorb significant energy from the electromagnetic field at those frequencies.

Figures 2(a) and 2(b) show the frequency dependence of the real and imaginary parts of the permittivity parallel (ε ) and perpendicular (ε ⊥ ) to the optical axis of α-quartz, respectively. The sharp resonance peaks of the imaginary parts of both ε and ε ⊥ , along with their small values for broad frequency ranges, show that the SPPs can propagate over relatively large distances within these frequency intervals. ε has a less number of absorbing peaks than ε ⊥ , and hence the SPP energy transport is expected to be more efficient along the optical axis (ε x = ε ) than along the perpendicular axis (ε x = ε ⊥ ).

The dispersion relation and propagation length of SPPs traveling along the film/air interface are shown in Figs. 3(a of SPPs. For SPPs traveling along the optical axis of α-quartz (ε x = ε and ε y = ε z = ε ⊥ ), this range is 101 Trad/s < ω < 233 Trad/s; while for SPP propagating along the axis of ε ⊥ (ε x = ε ⊥ and ε y = ε z = ε ), the frequency interval is 108 Trad/s < ω < 233 Trad/s. Note that for both propagation directions, β R decreases and approaches the light line, as the frequency decreases. This indicates a photon-like nature of the SPPs. As the frequency increases, the dispersion relation separates from the light line and tends to a phonon-like behavior. This feature strengths when the film thickness increases. The corresponding propagation lengths in Fig. 3(b) take larger values at the frequencies where the absorption of energy is minimal (Fig. 2). For each propagation direction, the maxima of Λ are independent of d and depends on both the real and imaginary parts of the respective permittivity, as was shown previously [6]. This is why the minimum of the propagation lengths are not placed exactly at the resonance frequency of ε I , but rather in its vicinity. Note also that for a given film thickness, larger propagation lengths are obtained when the optical axis is parallel to the SPP direction.

The real parts of the transverse wave vectors inside and outside of the thin film are positive, as shown in Fig. 4(a). This fact along with Figs. 3(b) reveals that the propagation of SPPs is present in a broad band of frequencies. This is further confirmed by the distribution of the Poynting vector shown in Fig. 4(b). The energy flux propagating inside the film is close to zero and negligible in comparison with that outside of it, due to the very small absorption by the surrounding air. The high concentration of energy at the interface enables the absorbing film to support the propagation of SPPs. As a result of the coupling of the electromagnetic fields traveling along both film interfaces, the energy flux increases as the film thickness decreases. This indicates that the the SPP thermal conductivity should increase as the film thickness reduces, which is consistent with Eq. ( 4) and Fig. 3 The integrand of the thermal conductivity defined in Eq. 4 as a function of the frequency is shown in Fig. 5(a), for SPP propagation along two principal axes of the αquartz film. Note that the oscillatory behavior of these spectra is determined by the permittivity, and their maxima coincides with the corresponding ones of the propagation lengths shown in Fig. 3(b). The major contributions to the SPP thermal conductivity (area under the curve) arises from the vicinity of these peaks, as shown by the cumulative thermal conductivity in Fig. 5(b). Note that the increase of the thermal conductivity exhibits an inflection point at the frequency where the maxima of Λ occur. The SPP thermal conductivity κ increases as the film thickness decreases, and it can be much larger than its phonon counterpart for a film thickness smaller than 170 nm, as shown in Fig. 5(c). It is clear that κ becomes significant for nano thin films and negligible for micro thin films, as reported previously [4,6]. κ is higher when the SPPs propagate along the principal axis with smaller absorption (the optical axis, in this case), as expected. The enhancement of κ, due to the permittivity anisotropy of α-quartz, increases for smaller film thicknesses and lower temperatures. For d = 100 nm, it passes from 25% at T = 300 K to 14% at T = 500 K. Thus, it is clear that κ is fairly sensitive to the propagation direction and temperature of the anisotropic film. In contrast to the phonon thermal conductivity of α-quartz, the SPP one increases with temperature (Fig. 5(d)). This shows that the SPPs have the potential to diminish the overheating of materials under high temperatures of operation, by increasing their thermal conductivity.

The source of the high SPP thermal conductivity is the propagation length at small wavevectors, which are very close to the light line, as shown in Figs. 3(a) and 3(b). The largest propagation length in this photon-like regime is about 1 m for a 100 nm-thick film of α-quartz. This long value of Λ implies that only thin films with a length (L) longer than this distance (L >> Λ) can benefit from the enhanced SPP thermal conductivity. In practice, however, L is not so long and therefore the propagation length has to be redefined. For this case, the Matthiessen's rule [START_REF] Ziman | Electrons and phonons: the theory of transport phenomena in solids[END_REF], establishes that the effective propagation length (Λ e ) is given by Λ -1 e = Λ -1 + L -1 . Thus, the effect of the film length L on the SPP energy transport appears through the reduction of the intrinsic propagation length Λ. To quantify the correspondin reduction of the SPP thermal conductivity, we use the explicit expression for κ in terms of the propagation length [4,6].

Figure 6(a) shows the effective propagation length of SPPs propagating along the optical axis of the thin film, for three lengths L. The case L > 2 m corresponds to Λ e = Λ. Note that Λ e follows the same trend than Λ and its maximum values are smaller but comparable to L. This reduction of the propagation length from meters to centimeters has a strong impact on the SPP thermal conductivity κ, as shown in Fig. 6(b). For T = 300 K, κ reduces from 82 W/m.K to 13 W/m.K when the film length is scaled down from L > 2 m to L = 10 cm. Despite of this sizeable reduction of κ, its value is only 1.1 W/m.K smaller than its phonon counterpart, which shows that the SPPs are still able to significantly enhance the phonon thermal conductivity of nano thin films. Even though lower values of κ are obtained for a smaller film length (L = 1 cm), this reduction can be offset at higher temperature, as shown in Fig. 6(c). In summary, the energy transport of SPPs along nano thin films can be optimized by choosing the propagation direction along the axis of less absorption, decreasing the film thickness, increasing the film length and raising the temperature. A simple expression for the SPP thermal conductivity as a linear combination of the inverse first and third powers of the film thickness has been derived. The thermal conductivity along the optical axis of a 100 nm-thick film of α-quartz is 25% higher than that along a perpendicular direction. A thermal conductivity of 13 W/m.K has been obtained for a thin film of 100 nm in thickness and 10 cm in length, at room temperature.

1 FIG. 1 .

 11 FIG.1. Cross-section of the thin film under consideration.

FIG. 2 .

 2 FIG. 2. Real (εR) and imaginary (εI ) parts of the permittivity (a) parallel and (b) perpendicular to the optical axis of αquartz as a function of frequency.

6 zFIG. 3 .

 63 FIG. 3. (a) Dispersion relation and (b) propagation length as a function of the frequency, for a thin film of α-quartz surrounded by air.

  (b). For the frequency ω = 169.9 Trad/s and thickness d = 200 nm, Figs. 3(b) and 4(b) show that the SPPs propagate about 100 mm along the thin film and just 100 µm across it. This difference on three orders of magnitude shows the promising potential of SPPs to transport energy along thin films.

FIG. 4 .

 4 FIG. 4. Frequency dependence of Re(pj) inside and outside of a 100 nm-thick film of α-quartz surrounded by air. (b) Poynting vector distribution at x = 0 and ω = 169.9 Trad/s, for SPPs propagating along the optical axis.

3 F

 3 FIG. 5. (a) Thermal conductivity spectrum and (b) its cumulative values for d = 100 nm and T = 300 K. Thermal conductivity as a function of the film (c) thickness and (d) temperature (d = 100 nm). The experimental points stand for the phonon counterpart of α-quartz along its optical axis.