Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Emmanuel Lorin
email: elorin@math.carleton.ca

Towards Perfectly Matched Layers for time-dependent space fractional PDEs

Keywords: fractional partial differential equations, unbounded domain, perfectly matched layer, Fourier pseudospectral approximation, time splitting scheme . . . . . . 39

Perfectly Matched Layers (PML) are proposed for time-dependent space fractional PDEs. Within this approach, widely used powerful Fourier solvers based on FFTs can be adapted without much effort to compute Initial Boundary Value Problems (IBVP) for well-posed fractional equations with absorbing boundary layers. We analyze mathematically the method and propose some illustrating numerical experiments.

Introduction

Over the last decade, the study of Fractional Partial Differential Equations (FPDEs) became a huge domain of investigation in many areas of science and engineering [START_REF] Guo | Fractional Partial Differential Equations and Their Numerical Solutions[END_REF][START_REF] Karniadakis | Special Issue on Fractional PDEs: Theory, Numerics, and Applications[END_REF][START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF]. Therefore, simulating accurately and efficiently a FPDE is a major challenge in computational science. In this paper, we are interested in suitably truncating the infinite spatial computational domain for one-dimensional (and higher-dimensional) space fractional Initial Value Problems (IVPs) of the form

∂ t u(t, x) + α∈R v [α] (x)∂ α x u(t, x) = f (u(t, x)), (1)
or

∂ t u(t, x) + α∈R v [α] (x)∆ α x u(t, x) = f (u(t, x)), (2)
where x ∈ R and t > 0, {v [α] } α∈R are smooth real-or purely complex-valued functions, R being a finite set of strictly positive real numbers, and f is square integrable. To this system is added an initial data u(t = 0, x) = u 0 (x). Spatially bounding the above IVPs naturally leads to Initial Boundary Value Problems (IBVPs) set on the truncated domain. Simple boundary conditions like homogeneous Dirichlet, Neumann or Robin boundary conditions as well as periodic boundary conditions can be used in some situations to get a finite computational domain. However, this choice of boundary condition does not always provide an admissible approximate solution. For integer order Partial Differential Equations (PDEs) (α ∈ N), many approaches can be used to build [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF] absorbing, transparent or non-reflecting boundary conditions (ABC, TBC, NRBC), according e.g. to the structure and physical/mathematical properties of the PDEs, the kind of expected scheme to be used... Usually, constructing such boundary conditions requires nontrivial mathematical analysis and can often be recast as the fundamental problem of building the exact Dirichlet-to-Neumann operator or some of its approximations for the associated PDE. In addition, several complex mathematical (definition, well-posedness) and computational issues (stability, accuracy, computational complexity) are usually related to such boundary conditions. For the space FPDEs, to the best of authors knowledge, only a few derivations of ABCs/TBCs/NRBCs have been recently obtained [START_REF] Baeumer | Boundary conditions for fractional diffusion[END_REF][START_REF] Cusimano | On reflecting boundary conditions for space-fractional equations on a finite interval: proof of the matrix transfer technique[END_REF][START_REF] Huang | Finite difference methods for fractional Laplacians[END_REF][START_REF] Kelly | Boundary conditions for two-sided fractional diffusion[END_REF]. Furthermore, closely related are newly designed truncation techniques for nonlocal models (see e.g. [START_REF] Du | Numerical solution of a two-dimensional nonlocal wave equation on unbounded domains[END_REF][START_REF] Du | Nonlocal wave propagation in unbounded multi-scale mediums[END_REF]), generalizing in some sense the notion of FPDEs. Let us remark that such boundary conditions were also nontrivially derived for time-fractional PDEs [START_REF] Li | An efficient second-order finite difference method for the one-dimensional Schrödinger equation with absorbing boundary conditions[END_REF][START_REF] Zhang | Efficient numerical computation of time-fractional nonlinear Schrödinger equations in unbounded domain[END_REF][START_REF] Zheng | Numerical solution of nonlocal diffusion equation on the real line[END_REF].

Finally, let us note that these boundary conditions cannot be a priori implemented into a Fourier pseudospectral approximation scheme [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Boyd | Chebychev and Fourier Spectral Methods[END_REF][START_REF] Li | Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[END_REF][START_REF] Shen | Spectral Methods: Algorithms, Analysis and Applications[END_REF] based on FFTs since the enforced boundary conditions are not periodic.

An alternative method to avoid spurious unphysical reflection at the domain boundary consists in using the method of Perfectly Matched Layers (PMLs) which has been extensively studied for many integer order PDE models (see e.g. [START_REF] Bécache | Perfectly mathched layers for the convected Helmholtz equation[END_REF][START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bérenger | Perfectly matched layer for the FDTD solution of wave-structure interaction problems[END_REF][START_REF] Bermúdez | An exact bounded perfectly matched layer for time-harmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF][START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF][START_REF] Hu | On absorbing boundary conditions for linearized Euler equations by a perfectly matched layer[END_REF][START_REF] Hu | A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF][START_REF] Turkel | Absorbing PML boundary layers for wave-like equations[END_REF][START_REF] Zeng | The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media[END_REF]). PMLs usually enjoy some good mathematical and numerical properties: easiness of implementation, flexibility, accuracy, stability,... which made them extremely popular since their introduction in the seminal paper by Bérenger [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] in 1994 for simulating electromagnetic waves. Even if PMLs are very attractive, to the best of our knowledge, the extension to the space FPDEs has never been studied. The aim of this paper is to contribute to deriving PMLs for some classes of FPDEs, resulting in Fractional PMLs (FPMLs). As for the integer case, these FPMLs are easy to integrate into the FPDE mathematical formulation. Of course, they also need to be fixed for each FPDE, as for the standard case, and in particular concerning the absorption profiles and their tuning parameters. Therefore, in the present paper, instead of focusing on one specific FPDE, we explain the general idea and illustrate the accuracy of FPMLs through explicit examples. For a better understanding of the full approach, more focused future investigations are required to develop optimal FPMLs for a given FPDE. In addition, we also introduce a specific Fourier based pseudospectral discretization scheme for approximating the FPMLs model which appears to be very flexible and accurate for the resulting modified problem. However, it is also clear that many other discretization approaches could be investigated, which should be again further studied. In particular, one can freely choose the boundary condition imposed at the outer domain boundary. Therefore, this allows for adapting many already existing discretization schemes. Even if we start with Riemann-Liouville operators [START_REF] Guo | Fractional Partial Differential Equations and Their Numerical Solutions[END_REF][START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF], the extension to other definitions of fractional derivative operators should be possible, as we notice for Caputo-type derivative operators. Moreover, we also prospect the extension to higher-dimensional problems and FPDEs involving fractional laplacians [START_REF] Kwasnicki | Ten equivalent definitions of the fractional laplace operator[END_REF][START_REF] Lischke | What is the fractional laplacian?[END_REF]. The method can be adapted but some nontrivial questions remain open for further improvements. Finally, well-posedness of the truncated problem is partially addressed here as well as a stability analysis of the scheme. Finally, extending the method to time-fractional PDE seems possible.

The paper is organized as follows. In Section 2, we present the derivation and possible approximations of PMLs for one-dimensional space FPDEs. The case of fractional laplacians is also considered. Section 3 is dedicated to the numerical approximation of the FPDE with FPML by a well-designed pseudospectral method combined with time-splitting schemes. Several numerical experiments illustrate the ideas and methods presented in this paper. Extension to higher dimensional problems is discussed in Section 4. In Section 5, we analyze the well-posedness of the FPDE with FPML, as well as the stability of an implicit scheme derived in this paper. We finally conclude in Section 6.

Fractional Perfectly Matched Layers (FPMLs) and their approximations

Some general definitions and basic properties about fractional derivatives are recalled below [START_REF] Di Nezza | Hitchhiker's guide to the fractional Sobolev spaces[END_REF][START_REF] Guo | Fractional Partial Differential Equations and Their Numerical Solutions[END_REF][START_REF] Karniadakis | Special Issue on Fractional PDEs: Theory, Numerics, and Applications[END_REF][START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF]. All along the paper, in particular in Section 5, some key results about fractional Sobolev spaces will be needed and recalled. Here, we simply provide the minimal requirement about fractional derivatives for a good understanding of the material presented in the paper.

Fractional Perfectly Matched Layers (FPMLs)

For α ∈ R * + and p = α + 1 ( α being the integer part of α), we define the special Gamma function Γ, for Re(z) > 0, as

Γ(z) := ∞ 0 x z-1 e x dx,
with Γ(p) = (p -1)!, for p ∈ N * . The fractional derivative that we consider in this paper is the so-called Riemann-Liouville (RL) fractional derivative of order α over the interval (-∞; x], which is defined by

RL D α x u(x) = 1 Γ(p -α) d p dx p x -∞ u(y) (x -y) α-p+1 dy, (3)
while the right RL fractional derivative is given by

RL x D α u(x) = (-1) p Γ(p -α) d p dx p +∞ x u(y) (y -x) α-p+1 dy .
Similarly, we introduce the left fractional Riemann-Liouville integral operator of order α as

RL I α x u(x) = 1 Γ(α) x -∞ u(y) (x -y) 1-α dy. (4)
Through these notations, we have the relation

RL D α x u(x) = d p dx p RL I p-α x u(x). (5)
Let us remark that many other fractional derivatives/integrals can also be defined [START_REF] Guo | Fractional Partial Differential Equations and Their Numerical Solutions[END_REF][START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF], like e.g. the well-known Caputo (C) derivative of order α

C D α x u(x) = 1 Γ(p -α) x -∞ u (p) (y) (x -y) α-p+1 dy. (6)
In the sequel of the paper, we restrict our study to the Riemann-Liouville fractional partial derivatives with respect to x, but it may be possible that some of the ideas can be extended to other definitions of the fractional operators as suggested below. From now on, the RL fractional derivative will be denoted by ∂ α x or x ∂ α , respectively corresponding to RL D α x and RL x D α . We again refer to [START_REF] Guo | Fractional Partial Differential Equations and Their Numerical Solutions[END_REF][START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF] for the mathematical properties of fractional operators.

Let us now consider the basic approach related to PMLs. The equations under consideration are defined on an open one-dimensional bounded physical domain denoted by D Phy , but the extension to higher-dimensional problems can be derived (see also Section 4). As usual, we add a layer which is called D PML , surrounding D Phys , stretching the x-coordinate. The overall computational domain is then defined by: D = D Phy ∪ D PML . For the onedimensional case, D = [-L, L] and D Phys = [-L * , L * ], with L * < L, the layer thickness being δ x . Usual PMLs require a complex stretching of the real spatial coordinate x such as

x(x) = x + e iθ x L * σ(s)ds, (7)
where the absorbing function σ : D → R is defined [START_REF] Antoine | Coupling spectral methods and perfectly matched layer for simulating the dynamics of nonlinear Schrödinger equations[END_REF] as

(α ∈ N * ) σ(x) = σ(|x| -L), L * |x| < L, 0, |x| < L * . (8)
The rotation angle θ is usually fixed by the problem under study. For example, θ = π/2 is often considered for (integer order) time harmonic Helmholtz-type problems [START_REF] Bécache | Perfectly mathched layers for the convected Helmholtz equation[END_REF][START_REF] Bermúdez | An exact bounded perfectly matched layer for time-harmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF][START_REF] Turkel | Absorbing PML boundary layers for wave-like equations[END_REF] while θ = π/4 is more adapted to Schrödinger problems [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Zeng | The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media[END_REF][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF].

Here, we derive the PML equation only for the right layer, the extension to the left side being straightforward by a symmetry argument. In addition, the boundary condition at point {L} can be fixed for example as a homogeneous Dirichlet or Neumann boundary condition or a periodic boundary condition, since the solution is supposed to be small enough after being damped. This provides a great flexibility for the choice of the discretization method. We will later choose a periodic boundary condition so that Fourier pseudospectral approximation methods based on FFTs can be used.

The starting point of PMLs for PDEs with integer order partial derivatives considers the following damped operator

∂ x → ∂ x := 1 (1 + e iθ σ(x)) ∂ x , (9)
modifying hence the initial PDE. Following this idea, if e.g. 0 < α < 1, we can introduce the new complex-valued fractional derivative operators as

∂ α x u(x) = 1 Γ(1 -α) 1 (1 + e iθ σ(x)) d dx x -∞ u(y) (x -y) α (1 + e iθ σ(y)) 1-α dy = 1 (1 + e iθ σ) ∂ α x ((1 + e iθ σ) 1-α u). (10)
For 1 < α < 2 (p = 2), then one gets from (3) and by a similar argument

∂ α x u(x) = 1 Γ(2 -α) 1 (1 + e iθ σ(x)) d dx ( 1 (1 + e iθ σ(x)) d dx x -∞ u(y) (x -y) α-1 (1 + e iθ σ(y)) 2-α dy) = 1 (1 + e iθ σ) d dx ( 1 (1 + e iθ σ) d dx RL I 2-α x ((1 + e iθ σ) 2-α u))), (11)
which degenerates, by evaluating from the inner to the outer parentheses, into the standard formula

∂ 1 xu(x) = 1 Γ(1)
d 2 dx 2 x -∞ u(y)dy = 1 1 + e iθ σ du dx , (12)
that could be also obtained from [START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF] taking the lower limit on α. Similarly, one gets

∂ 2 xu = 1 Γ(2) d 2 dx 2 u = 1 (1 + e iθ σ) d dx ( 1 (1 + e iθ σ) du dx ). ( 13
)
A generalization of the absorbing RL fractional derivative operators of order α ∈ R * could be derived following the formula

∂ α x u(x) = d p dx p RL I p-α x ((1 + e iθ σ) p-α u), (14)
setting

d p dx p = 1 (1 + e iθ σ) d dx ( 1 (1 + e iθ σ) d dx (...)) (p times).
In the case of the Caputo derivative [START_REF] Baeumer | Boundary conditions for fractional diffusion[END_REF], then a natural extension of a PML Caputo fractional derivative operator of order α ∈ R * is

C D α x u(x) = 1 Γ(p -α) x -∞ u (p) (y) (x -y) α-p+1 (1 + e iθ σ) -α dy = RL I p-α x ((1 + e iθ σ) -α u (p) ). (15)
Choosing an absorption function is not trivial and is based on the underlying PDE model. The first profiles proposed in the literature are based on polynomial functions. Here, we consider the quadratic and cubic functions Type I:

σ 0 (x + δ x ) 2 , Type II: σ 0 (x + δ x ) 3 , ( 16
)
where σ 0 is a real-valued positive parameter to adjust. Other possibilities derived for example for the Helmholtz equation [START_REF] Bermúdez | An exact bounded perfectly matched layer for time-harmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF] and later used for the Schrödinger equation [START_REF] Antoine | Coupling spectral methods and perfectly matched layer for simulating the dynamics of nonlinear Schrödinger equations[END_REF] include unbounded functions σ such that

L L * σ(x)dx = +∞. ( 17
)
Among these, we use the following discontinuous and continuous functions Type III: -σ 0 x , Type IV:

σ 0 x 2 , Type V: - σ 0 x - σ 0 δ x , Type VI: σ 0 x 2 - σ 0 δ 2 x . (18)
Let us remark that other choices of absorbing functions may be more suited to fractional equations, but this point is not analyzed in the current paper. It is well-known that tuning the parameter σ 0 in PMLs is relatively problem dependent. As noticed e.g. in [START_REF] Antoine | Coupling spectral methods and perfectly matched layer for simulating the dynamics of nonlinear Schrödinger equations[END_REF][START_REF] Bermúdez | An exact bounded perfectly matched layer for time-harmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF], Types III-VI usually exhibit a larger stability region for a numerical discretization for Helmholtz and Schrödinger problems.

Approximate FPMLs for RL operators

Even if the above Fractional PML (FPML) operators could be used, we propose here an approximate and modified version of a PML fractional derivative operator which can be related to the previous RL operators. Let us first recall that, for the RL operators, we have the following Fourier symbols

F(∂ α x u)(ξ) = (iξ) α F(u)(ξ), F( RL I α x u)(ξ) = (iξ) -α F(u)(ξ), (19)
where F is the Fourier transform according to x and ξ is the corresponding dual co-variable. Therefore, if one considers the principal symbol σ p in the sense of pseudodifferential operators [START_REF] Taylor | Pseudodifferential Operators[END_REF] of the operator ∂ α x defined by [START_REF] Boyd | Chebychev and Fourier Spectral Methods[END_REF], one gets, for α ∈ R * ,

σ p (∂ α x ) = (1 + e iθ σ) -α (iξ) α . (20)
Consequently, another fractional PML operator that could be used to bound the computational domain is

∂ α x := 1 (1 + e iθ σ) α ∂ α x , (21)
then generalizing the well-know formula [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] for integer order derivatives to fractional orders of differentiation. Different strategies can be developed for constructing the absorbing operators. Even if ( 21) is a reasonable approximation for α < 1, this is not the unique choice.

For α > 1, we would like to mimic this for the FPML [START_REF] Boyd | Chebychev and Fourier Spectral Methods[END_REF] corresponding to α = p, whenever p -1 < α < p. Within this aim, we set, for any α ∈ R * + ,

α = n α + q α ∈ R * + , β α := α n α + 1 ∈ R * + , (22)
with

n α := α -1, if α ∈ N, α , if α ∈ R -N, ( 23
)
where α is the integer part of α, and q α := α -n α ∈ R * + is strictly less than 1. The main idea consists in decomposing the derivatives ∂ α

x , with α ∈ R * + , as

σ(∂ α x ) = (iξ) α = ((iξ) βα ) nα+1 = (iξ) qα (iξ) nα , ∂ α x = (∂ βα x ) nα+1 = ∂ qα x ∂ nα x .
We then consider the following various FPML operators (which are approximations up to some lower order operators of the original FPML operators), for α > 1, where β α ∈ (0, 1), for ∂ α x

• Approach 0. For n α defined by [START_REF] Du | Nonlocal wave propagation in unbounded multi-scale mediums[END_REF], the first approach corresponds to (11)

∂ α x u = 1 (1 + e iθ σ) d dx ( 1 (1 + e iθ σ) • • • d dx RL I nα-α x ((1 + e iθ σ) nα+1-α u)).
• Approach 1. This approach directly mimics [START_REF] Boyd | Chebychev and Fourier Spectral Methods[END_REF], then leading to

σ(∂ α x ) ≈ (1 + e iθ σ) -βα (iξ) βα nα+1 .
For u smooth enough, this leads to the corresponding approximation

∂ α x u ≈ ( 1 (1 + e iθ σ βα ∂ βα x ) nα+1 u. (24)
• Approach 2a. A closely related approach consists in approximating the symbol σ(∂ α x ) as follows

σ(∂ α x ) ≈ (1 + e iθ σ) -1 (iξ) qα (1 + e iθ σ) -1 (iξ) nα , yielding ∂ α x u ≈ 1 1 + e iθ σ qα ∂ qα x ( 1 1 + e iθ σ∂ x nα u. (25)
• Approach 2b. Symmetrically, one can consider

∂ α x u ≈ 1 (1 + e iθ σ) nα ∂ nα x ( 1 1 + e iθ σ qα ∂ qα x )u. ( 26
)
• Approach 3. As previously explained, the last approach proposes to select the principal symbol to ∂ α x , i.e.

∂ α x u ≈ 1 1 + e iθ σ α ∂ α x u. (27)
In Approaches 1 and 2, if α = p ∈ N -{0, 1}, the RHS of ( 24) and ( 25) becomes [START_REF] Boyd | Chebychev and Fourier Spectral Methods[END_REF].

Example 1. As an illustration and to compare the various strategies, we report |u| (in logscale) at T = 10 on Fig. 1 for the different FPML approaches, when solving the FPDE

∂ t u(t, x) -v(x)∂ 3/2 x u(t, x) = 0, t > 0, x ∈ R, u(t, 0) = u 0 (x), x ∈ R, (28)
where u 0 (x) = N -1 e -2x 2 +ik 0 x (N being a L 2 (D)-normalization constant), for k 0 = 1, the multiplicative function being v(x) = e -5×10 -2 x 2 . The physical domain of computation is D Phys := (-L * ; L * ) while the extended domain with FPMLs is D := (-L; L), setting L = 8 and L * = 0.95L. The approximate FPMLs is chosen using the Type I absorbing function [START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF], with σ 0 = 10, θ = π/8, and periodic boundary conditions at the endpoints of the layers.

More specifically, we are solving the IBVP

   ∂ t u(t, x) -v(x)∂ 3/2 x u(t, x) = 0, t > 0, x ∈ D, u(t, -L) = u(t, L), t > 0, u(t, 0) = u 0 (x), x ∈ D, (29)
where ∂ 3/2 x is one of the FPML operators. Since α = 3/2, then we have: n α = 1, q α = 1/2 and β α = 3/4. The reference solution is computed on a larger domain so that no reflection arises at the right boundary. In addition, we also report the solution computed on D but with periodic boundary conditions and no FPML. The discretization parameters are: N x = 1001 and ∆t = 10 -2 (N t = 2500). The numerical method is described later in the paper (Section 3). Let us remark that any other discretization could be used for the FPML method, the scheme depending on the boundary conditions set at the endpoints. This first example shows that the different FPML approaches have relatively similar absorption features, which are all very effective. In particular, and to fix the ideas in the paper, we consider now the FPML based on the operator proposed in Approach 1, except if specified otherwise, which is well-suited in the numerical simulation while being simple. Since numerous parameters need to be tuned for optimizing the FPML (e.g. Type of PML, σ 0 , θ, δ x ) and because this is problem dependent, we will address in details the question of selecting the best absorbing operator approximation for specific FPDE in a future work.

Focus on FPMLs for the fractional laplacian

Let us now specifically look at the problem of building FPMLs related to FPDEs involving the 1D fractional laplacian α x , 0 < α 1, which appears in many fractional models [START_REF] Baeumer | Boundary conditions for fractional diffusion[END_REF][START_REF] Du | Numerical solution of a two-dimensional nonlocal wave equation on unbounded domains[END_REF][START_REF] Du | Nonlocal wave propagation in unbounded multi-scale mediums[END_REF][START_REF] Guo | Fractional Partial Differential Equations and Their Numerical Solutions[END_REF][START_REF] Karniadakis | Special Issue on Fractional PDEs: Theory, Numerics, and Applications[END_REF][START_REF] Lischke | What is the fractional laplacian?[END_REF]. Let us remark here that we do not use the operator (x) α and rather include the normalization constant e iπα in the multiplicative function v if necessary. First, we notice that For u smooth enough, this leads to the following possible choices of operators to define the FPMLs Approach 0.

σ( α x ) = e iαπ |ξ| 2α = e iβαπ |ξ| 2βα ) nα+1 = e iqαπ |ξ| 2qα × e inαπ |ξ| 2nα , α x = ( βα x ) nα+1 = qα x nα x .
α x u = 1 (1 + e iθ σ) 1/2 x • • • 1 (1 + e iθ σ) 1/2 x RL I nα-α x ((1 + e iθ σ) nα+1-α u), Approach 1. α x u ≈ 1 (1 + e iθ σ βα βα/2 x 2(nα+1) u, Approach 2a. α x u ≈ 1 (1 + e iθ σ) qα qα/2 x 1 (1 + e iθ σ 1/2 1/2 x 2(nα+1) u, Approach 2b. α x u ≈ 1 (1 + e iθ σ 1/2 1/2 x 2(nα+1) 1 (1 + e iθ σ) qα qα/2 x u, Approach 3. α x u ≈ 1 1 + e iθ σ 2α α x u.
(30) Notice that, for the Approach 1 and for ∂ 2α

x , one considers

σ(∂ 2α x ) ≈ (1 + e iθ σ) -β 2α (iξ) β 2α n 2α +1 . For α in R -N * , we have β 2α = 2α/(1 + 2α ) and n 2α = 2α . Then, one gets the approximation σ(∂ 2α x ) ≈ (1 + e iθ σ) -2α 1+ 2α (iξ) 2α 1+ 2α
2α +1 , leading to the approximate FPML operator

∂ 2α x u ≈ 1 (1 + e iθ σ) 2α 1+ 2α ∂ 2α 1+ 2α x 2α +1 u = 1 (1 + e iθ σ) β 2α ∂ β 2α x n 2α +1 u. (31)
As an example, from ( 30) and ( 31) and taking α = 3/4, n α = 0,

β α = 3/4, n 2α = 1, β 2α = 3/4, one gets 3/4 x u ≈ 1 (1 + e iθ σ) 3/4 3/8 x 2 u, ∂ 3/2 x u ≈ 1 (1 + e iθ σ) 3/4 ∂ 3/4 x 2 u.
Example 2. To compare the different FPMLs, we plot |u| (in logscale) at T = 10 on Fig. 2 for several choices of FPMLs when solving the FPDE (linear fractional Schrödinger-type equation)

i∂ t u(t, x) + v(x)∆ 3/4 x u(t, x) = 0, t > 0, x ∈ R, u(t, 0) = u 0 (x), x ∈ R, (32)
with

u 0 (x) = N -1 e -2x 2 +ik 0 x (where N is the L 2 (D)-normalization constant), the function v is v(x) = e -5×10 -2 x 2 +3iπ/4
and the wave number is k 0 = 5. The physical domain of computation is D Phys := [-L * ; L * ] while the entire domain with FPMLs is D := (-L; L), choosing L = 8 and L * = 0.95L. Regarding the FPMLs, we fix the Type I absorbing function ( 16), for σ 0 = 10, θ = π/8, and we consider periodic boundary conditions, i.e., we solve for example the IBVP

       i∂ t u(t, x) + v(x) 1 (1 + e iθ σ(x)) 3/4 3/8 x 2 u(t, x) = 0, t > 0, x ∈ D, u(t, -L) = u(t, L), t > 0, u(t, 0) = u 0 (x), x ∈ D. (33)
The discretization parameters are again N x = 1001 and ∆t = 10 -2 (N t = 1000). We conclude that all the FPMLs are highly absorbing and relatively similar in terms of quality (up to the optimization of the tuning parameters). From now on, we fix the FPML based on Approach 1.

FPMLs pseudospectral approximation schemes: implementation and numerical examples

In this section, we propose a general methodology for solving IBVPs using Fourier-based pseudospectral methods [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Boyd | Chebychev and Fourier Spectral Methods[END_REF][START_REF] Li | Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[END_REF][START_REF] Shen | Spectral Methods: Algorithms, Analysis and Applications[END_REF] with FPML. Naturally, other approximation schemes could also be used. We start with fractional equations involving spatial orders of derivation r < 1 in Subsection 3.1, and next consider the more general case r ∈ R * + in Subsection 3.2. Initially, the problem of interest is an IVP set in R. For obvious computational reasons, the problem must be solved as an IBVP on a truncated domain, as it is commonly the case. Considering the well-posedness, suitable boundary conditions need to be imposed at the boundary of the truncated domain. When Fourier-based methods are used, periodic boundary conditions are naturally fixed, which are then potentially inappropriate for computing delocalized solutions. The proposed computational methodology developed here allows for fixing this issue in the case of FPDEs, based on the previous FPMLs combined with a pseudospectral discretization and time splitting schemes.

IBVP with space fractional operators of order α less than 1

We consider the following IBVP, which is assumed to be well-posed in [0; T ] × D Phy (see Section 5), where the equation contains fractional operators of order α strictly less than 1.

For a given initial data u(t = 0, •) = u 0 that we assume to be compactly supported in D Phy , we want to solve the IVP

∂ t u(t, x) + α∈R v [α] (x)∂ α x u(t, x) = f (u(t, x)),
where {v [α] } α∈R are smooth functions, f is square integrable, and R is a finite set of strictly positive real numbers such that the IBVP is well-posed, for 0 < r := max R < 1. Generalizing the idea presented in Subsection 2.2 for the Approach 1, we set

S(x) = 1, if |x| < L * , 1 + e iθ σ(x), if L * |x| < L, (34)
and consider the modified IBVP on [0; T ] × D with FPMLs

∂ t u(t, x) + α∈R v [α] (x) S α (x) ∂ α x u(t, x) = f (u(t, x)), (35)
with periodic boundary conditions u(t, -L) = u(t, L) and initial data u(t = 0, •) = u 0 .

We will also study the case of FPDEs of the form

∂ t u(t, x) + α∈R v [α] (x) α x u(t, x) = f (u(t, x)),
which leads to the new system

∂ t u(t, x) + α∈R v [α] (x) S α (x) α x u(t, x) = f (u(t, x)),
with initial data u(t = 0, •) = u 0 , and periodic boundary conditions at the domain interface {-L; L}. Let us denote the set of grid-points by

D Nx = x k 1 k 1 ∈O Nx , O Nx = k 1 ∈ N/ k 1 = 0, • • • , N x -1 ,
and the uniform mesh size by h = h

x := x k 1 +1 -x k 1 = 2L/N
x (for the entire domain D). The corresponding discrete wave numbers are defined by

ξ p = pπ/L, for p ∈ {-N x /2, • • • , N x /2- 1}.
Concerning the approximations, we use the following notation

u p (t) = Nx-1 k 1 =0 u k 1 (t)e -iξp(x k 1 +L) , u k 1 (t) = 1 N x Nx/2-1 p=-Nx/2 u p (t)e iξp(x k 1 +L) . (36)
Finally, we define the Fourier-based approximations of the fractional derivative operators of order α > 0 as

∂ α x u(t n , x k 1 ) ≈ [[∂ α x ]]u k 1 (t n ) := 1 N x Nx/2-1 p=-Nx/2 iξ p α u p (t n )e iξp(x k 1 +a) . (37)
Similarly, for the fractional laplacian, we introduce

α x u(t n , x k 1 ) ≈ [[ α x ]]u k 1 (t n ) := 1 N x Nx/2-1 p=-Nx/2 ξ p 2α e iαπ u p (t n )e iξp(x k 1 +a) . (38)
We denote by u n h the approximation of u at time t n on D Nx , and by S h and v

[α]
h the respective values of S and v [α] evaluated on the grid set D Nx . A simple second-order semi-explicit operator splitting in time of [START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF] reads

                 u n+1/2 h = u n h + ∆t 2 f u n+1/2 h , u n * h = u n+1/2 h -∆t α∈R v [α] h S α h [[∂ α x ]]u n+1/2 h , u n+1 h = u n * h + ∆t 2 f u n+1 h .
In principle, the above scheme is constrained by a CFL condition which limits its interest and application, in particular when the equation involves high-order (fractional) derivatives.

We conjecture that the above one-dimensional scheme is stable under a CFL condition of the form ∆t/∆x r C, for some C > 0 [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations[END_REF]. The larger r < 1, the more restrictive the above CFL condition. To avoid this problem, we rather consider an implicit version of the splitting scheme

                 u n+1/2 h = u n h + ∆t 2 f u n+1/2 h , u n * h = u n+1/2 h -∆t α∈R v [α] h S α h [[∂ α x ]]u n * h , u n+1 h = u n * h + ∆t 2 f u n+1 h . (39)
The scheme [START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF] requires the use of a linear system solver for the second fractional step, but ensures in principle better L 2 -stability of the scheme. Regarding the spatial discretization of the second equation in [START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF], from t n to t n * = t n + ∆t with initial data

u n+1/2 h := u n+1/2 j j∈O Nx
, it is equivalent to the linear system

A n h u n * h = u n+1/2 h , (40)
where

A n h = I + ∆t α∈R v [α] h S α h [[∂ α x ]], (41)
I being the identity matrix. From an implementation and computational point of view, this approach seems to be more complex, since the equation involves space variable coefficients. However, the system solution can be obtained thanks to the combination of a Krylov subspace iterative solver [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] like the GMRES with matrix-free evaluations based on FFTs. In addition, analytical preconditioners could be added to improve the convergence rate. We refer to [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF] for the introduction of this idea for the Gross-Pitaevskii equation in the framework of Bose-Einstein condensates which can be directly extended here to FPDEs. In the present paper, we use the GMRES without preconditioner [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF].

Example 3. To analyze how the proposed scheme is working, we consider the following linear FPDE

∂ t u -v [9/10] ∂ 9/10 x u = 0, t > 0, x ∈ R, u 0 (x) = N -1 e -(x-5) 2 +ik 0 x , x ∈ R,
with N a normalizing coefficient, v [9/10] (x) = e -x 2 /200 /2 and k 0 = -1. For the FPML, we use the Type I absorbing function given by [START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF], where σ 0 = 5 × 10 2 and θ = π/8, and with a layer thickness given by δ x := 0.05 for the domain D = [-10, 10]. The modified FPDE problem with FPML is then

       ∂ t u(t, x) - v [9/10] (x) S 9/10 (x) ∂ 9/10 x u(t, x) = 0, t > 0, x ∈ D, u(t, -L) = u(t, L), t > 0, u 0 (x) = N -1 e -(x-5) 2 +ik 0 x , x ∈ D, (42)
which is approximated by

I -∆t v 9/10 h S 9/10 h [[∂ 9/10 x ]] u n+1 h = u n h . ( 43
)
The discretization parameters are N x = 1001 and ∆t = 5 × 10 -2 . We plot in Fig. 3 the solutions {(x, t, log |u(x, t)|), (x, t) ∈ D × [0; T ]} for the FPML of Type I, for the solution with pure periodic boundary conditions (no FPML) and for a reference solution (computed on the larger domain [-20, 20]). We observe that the FPML properly works, and that there is almost no reflection at the boundary. A small amount of the solution passes through the boundary and comes back on the left boundary because of the periodic boundary condition. Example 4. Very similar ideas can be adapted to the case of equations involving fractional laplacians with FPMLs. Rather than presenting the general theory, we illustrate this through an example considering the following nonlinear FPDE

i∂ t u + v [1/2] 1/4 x u + v [3/4] 3/8 x u = f (u), t > 0, x ∈ R, u 0 (x) = N -1 e -(x+2) 2 +ik 0 x , x ∈ R,
where v [1/2] = e -x 2 /20+iπ/4 , v [3/4] (x) = e -x 2 /10+3iπ/8 , f : u → iκ|u| 2 u for κ = 10, and k 0 = 10. To include the FPML on D := [-L; L], the following modified problem is solved

∂ t u + v [1/2] S 1/2 1/4 x u + v [3/4] S 3/4 3/8 x u = f (u), (44)
with initial data and periodic boundary conditions. It is approximated as follows

               u n+1/2 h = u n h + iκ ∆t 2 u n h 2 u n+1/2 h , I + ∆t v [1/2] h S 1/2 h [[ 1/4 x ]] + ∆t v 3/4 h S 3/4 h [[ 3/8 x ]] u n * h = u n+1/2 h , u n+1 h = u n * h + iκ ∆t 2 u n * h 2 u n+1 h .
Alternatively, the ODEs involving the source term can be treated semi-analytically, i.e.

               u n+1/2 h = u n h exp i∆tκ ∆t 2 u n-1 h 2 , I + ∆t v [1/2] h S 1/2 h [[ 1/4 x ]] + ∆t v 3/4 h S 3/4 h [[ 3/8 x ]] u n * h = u n+1/2 h , u n+1 h = u n * h exp i∆tκ ∆t 2 u n * h 2 .
For the simulation, we fix ∆t = 5 × 10 -2 and N x = 1001 for the discretization. Regarding the FPML parameters, we choose L = 4, δ x = 0.075L and σ 0 = 10 for the Type I absorbing function. We compare {(x, t, log |u(x, t)|), (x, t) ∈ D × [0; T ]} in Fig. 4 the solution with periodic boundary conditions (no FPML), with FPML and for a reference solution (computed on the domain, [-16, 16]). This clearly illustrates the accuracy of the FPML combined with the splitting scheme and pseudospectral approximation. In the second part of the test, we compare the solution with three different types of FPML: Type I with σ 0 = 10, Type II with σ 0 = 10 2 , Type VI with σ 0 = 2.5 × 10 -2 (see Fig. 5). The FPML of Types I & II have an accurate absorbing effect, while Type VI-PML is less efficient for this example. We also compare on Fig. 6 the absorbing effect of the FPML of Type I with different values of σ 0 , i.e. σ 0 := 10 -1 , 1, 10, 10 2 , 10 3 . The best absorptions are obtained for the values σ 0 = 10 and 100.

IBVP with space fractional operators of order α larger or equal to 1

We assume now that the IBVP contains fractional operators of order larger or equal to 1. We impose u(t = 0, •) = u 0 , and consider

∂ t u(t, x) + α∈R v [α] (x)∂ α x u(t, x) = f (u(t, x)), (45)
where {v [α] } α∈R and f are smooth functions, and R is a finite set of positive real numbers such that the IBVP is well-posed. In the following, we propose to analyze the Approaches 1 and 2 from Subsection 2.2. Basically, the main idea consists in decomposing the derivatives ∂ α x , with α ∈ R * + . We can rewrite the equation ( 45) as

∂ t u(t, x) + α∈R v [α] (x)(∂ α nα+1 x ) nα+1 u(t, x) = f (u(t, x)) . ( 46
)
We then consider the modified IBVP on D including the FPML, following Approach 1 from Subsection 2.2,

∂ t u(t, x) + α∈R v [α] (x) 1 S βα (x) ∂ βα x nα+1 u(t, x) = f (u(t, x)), (47)
with periodic boundary conditions at ∂D and compactly supported initial data u 0 . By using the same notations as above, we have the following approximation where β α is defined in [START_REF] Du | Numerical solution of a two-dimensional nonlocal wave equation on unbounded domains[END_REF]. The implicit scheme reads

∂ βα x u(t n , x k 1 ) ≈ [[∂ βα x ]]u k 1 (t n ) := 1 N x Nx/2-1 p=-Nx/2 iξ p βα u p (t n )e iξp(x k 1 +L) , (48)
               u n+1/2 h = u n h + ∆t 2 f u n+1/2 h , u n * h = u n+1/2 h -∆t α∈R v [α] h 1 S βα h [[∂ βα x ]] nα+1 u n * h , u n+1 h = u n * h + ∆t 2 f u n+1 h . (49)
Approach 2a consists in solving the modified IBVP with FPML

∂ t u(t, x) + α∈R v [α] (x) 1 S qα (x) ∂ qα x 1 S(x) ∂ x nα u(t, x) = f (u(t, x)), (50)
where q α is defined in [START_REF] Du | Numerical solution of a two-dimensional nonlocal wave equation on unbounded domains[END_REF]. The corresponding second-order implicit splitting scheme is given by

                 u n+1/2 h = u n h + ∆t 2 f u n+1/2 h , u n * h = u n+1/2 h -∆t α∈R v [α] h S qα h [[∂ qα x ]] 1 S h [[∂ x ]] nα u n * h , u n+1 h = u n * h + ∆t 2 f u n+1 h . (51)
Practically, the second equation in ( 49) and ( 51) are linear systems solved by using the GMRES or BiCGStab iterative solvers [START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF] applied e.g. to

I + ∆t α∈R v [α] h S qα h [[∂ qα x ]] 1 S h [[∂ x ]] nα u n * h = u n+1/2 h .
We use here a GMRES algorithm to solve these linear systems.

Similarly, for the fractional laplacian and for the initial data u(t = 0, •) = u 0 , we consider

∂ t u(t, x) + α∈R v [α] (x) α x u(t, x) = f (u(t, x)), (52)
where {v [α] } α∈R and f are regular functions. The set R is a finite sequence of positive real numbers so that the IBVP is well-posed. Applying the Approach 1 from Subsection 2.2 yields

∂ t u(t, x) + α∈R v [α] (x) 1 S βα (x) βα/2 x nα+1 u(t, x) = f (u(t, x)),
while Approach 2 leads to

∂ t u(t, x) + α∈R v [α] (x) S(x) qα qα/2 x 1 S(x) 1/2 1/2 x 2(nα+1) u(t, x) = f (u(t, x)),
for t > 0 and x ∈ D. For well-posedness, periodic boundary conditions are added.

Example 5. In this example, the linear fractional system under consideration is

∂ t u(t, x) -v(x)∂ 3/2 x u(t, x) = f (u(t, x)), t > 0, x ∈ R, u(t, 0) = u 0 (x), x ∈ R, (53)
where the initial data is u 0 (x) = N -1 e -5×10 -1 (x+20/3) 2 +ik 0 x , with k 0 = -1. The multiplicative function is the gaussian v(x) = e -5×10 -3 x 2 and the nonlinear cubic term is f (u) = 10i|u| 2 u.

A FPML is next introduced into the formulation, with periodic boundary conditions. The computational domain is D = [-L; L], for L = 10 and L * = 0.95L. Setting α = 3/2, we have: n α = 1, q α = 1/2 and β α = 3/4. If one chooses the Approach 1 (see Eq. 47), we get the FPML equation

∂ t u(t, x) - v(x) S βα (x) ∂ βα x 1 S βα (x) ∂ βα x u(t, x) = f (u(t, x)),
while for Approach 2a (see Eq. 50) we have

∂ t u(t, x) - v(x) S qα (x) ∂ qα x 1 S(x) ∂ x u(t, x) = f (u(t, x)).
We fix ∆t = 5 × 10 -3 and N x = 2001 for the discretization.

In Fig. 7, we report |u(x, t)| (logscale) on D × [0; T ] (with T = 6), for respectively the periodic solution without FPML, the solution with FPML (Approach 2a) of Type I (with σ 0 = 10, θ = π/8), and a reference solution (computed on [-20, 20]). We clearly see that the FPMLs with pseudospectral approximation performs very well. For completeness, we also provide in Fig. 8 a zoom of the solution on D × [T /2; T ], where we see that the FPML solution disperses a little bit. We next compare in Fig. 9

i∂ t u = -( 1/2 x ) 2 u + κ|u| 2 u, t > 0, x ∈ R. ( 54
)
that we write as a FPDE (since x := (

1/ 2
x ) 2 ) to show that the proposed formulation can also work well even in this situation. The modified FPML equation [START_REF] Antoine | Coupling spectral methods and perfectly matched layer for simulating the dynamics of nonlinear Schrödinger equations[END_REF] reads with S given by Eqs. [START_REF] Lischke | What is the fractional laplacian?[END_REF]. For stability reasons, we still use the following implicit splitting scheme

i∂ t u = - 1 S 1/2 x 1 S 1/2 x u, (55)
             u n+1/2 h = u n h -i κ∆t 2 u n h 2 u n+1/2 h , u n * h = u n+1/2 h + i∆t 1 S h [[ 1/2 x ]] 1 S h [[ 1/2 x ]] u n * h , u n+1 h = u n * h -i κ∆t 2 u n * h 2 u n+1 h .
The initial data is the centered normalized Gaussian u 0 (x) = N -1 e -x 2 /2+ik 0 x , with wave number k 0 = 5. We fix κ = 10. The FPML of Type I, with σ 0 = 5 × 10 1 , θ = π/8 and δ x = 5 × 10 -2 L is used. The discretization parameters are fixed to ∆t = 10 -3 and N x = 1001 for the computational domain D = [-8, 8]. In Fig. 10, we report {(x, t, log |u(x, t)|), (x, t) ∈ D ×[0; T ]} for the periodic boundary conditions case (no FPML), for the FPML solution and the reference solution (computed on [-16, 16]). This illustrates the property that the FPML and its implementation is very efficient and accurate. We draw in Fig. 11 the amplitude |u| (logscale) of the three solutions at time T = 1.25. Let us remark here that the FPML that we use for the integer order PDE does not match with the standard way of writing a PML for Schrödinger-like PDEs (and Helmholtztype equations). Indeed, usually, the PML modification of the laplacian [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] is

∆ x → 1 S ∂ x 1 S ∂ x , (56)
while here, we have

∆ x → 1 S 1/2 x 1 S 1/2 x . ( 57
)
For the transformation (56), the involved operators are local and based on ∂ x while

1/2 x
appearing in (57) corresponds to a nonlocal operator. This operator is nontrivial to numerically approximate if we are not using a Fourier pseudospectral approximation scheme. The possibility of writing two kinds of PMLs for the integer case is related to the fact that the symbol of -∆ x is |ξ| 2 , which can be written as |ξ| × |ξ| (nonlocal) but also ξ × ξ (local).

Example 7. The equation under consideration is now the fractional nonlinear cubic Schrödinger equation

i∂ t u + v 3/4 x u = f (u), t > 0, x ∈ R,
with initial data u 0 (x) = N -1 e -5×10 -1 (x+4) 2 +ik 0 x (for k 0 = 5). The function v is given by v(x) = e -5×10 -2 x 2 +3iπ/4 and f (u) = 10i|u 2 |u. For the FPML, we fix the computational domain to D = [-L; L], with L = 8 and L * = 0.85L. We set α = 3/2, leading to n α = 1, q α = 1/2 and β α = 3/4. According to the methods developed above, we can either solve: Approach 1

i∂ t u + v S 3/4 3/8 x 1 S 3/4 3/8 x u = f (u), t > 0, x ∈ D, or Approach 2a i∂ t u + v S 1/2 1/4 x 1 S 1/2 x u = f (u), t > 0, x ∈ D,
with initial data u 0 and periodic boundary conditions. We fix ∆t = 10 -2 , N x = 501, and T = 10. In Fig. 12, we plot {(x, t, log |u(x, t)|), (x, t) ∈ D × [0; T ]} for the periodic boundary conditions (no PML), with the FPML of Type II (with σ 0 = 10, θ = π/8) and finally the reference solution (computed on [-16, 16]). This again shows that the FPMLs are very accurate when implemented in the Fourier pseudospectral method with time splitting. We also compare on Fig. 13 (Left) the absorption effects of the FPMLs with three different types of absorption function: Type I, Type II (with σ 0 = 10), and Type VI (with σ 0 = 5 × 10 -2 ), for δ x = 0.05L and θ = π/8. A similar test was performed where we compare the Approaches 1 and 2 on Fig. 13 (Right) with FPML of Types I & II (σ 0 = 10). The comparisons show a similar behavior for both approaches for the FPML with Type II absorption function, but the second approach seems more appropriate for the FPML with Type I function. Similar tests were also performed with an explicit scheme. However, the latter requires very small time steps for stability reasons and therefore implicit schemes are strongly recommended (see Section 5). Remark 3.1. The method developed in this paper can straightforwardly be extended to wellposed quasilinear equations

∂ t u + α∈R v [α] (u)∂ α x u = f (u), t > 0, x ∈ R.
For the sake of simplicity, let us assume that f = 0. For example, for Approach 1 in Subsection 2, we first rewrite the equation on D, including the FPML, as

∂ t u + α∈R v [α] (u) 1 S βα (x) ∂ βα x nα+1 u = 0 ,
where S is defined by [START_REF] Lischke | What is the fractional laplacian?[END_REF] and periodic boundary conditions are added. For the implicit scheme, we have

u n+1 h = u n h -∆t α∈R v [α] (u n h ) 1 S βα h [[∂ βα x ]] nα+1 u n+1 h .

Higher dimensional FPDEs

Approximate FPMLs for RL operators

We extend now the ideas and methods to two-dimensional problems (the 3D case can be treated by a direct extension). The numerical solution of the equation under consideration is approximated on an open two-dimensional bounded rectangular physical domain denoted by D Phy . As usual, we add a layer, denoted by D PML , surrounding D Phys , stretching the ν-coordinates, with ν = x, y. The overall computational domain is then:

D = D Phy ∪ D PML . In two dimensions, one gets D = [-L x , L x ] × [-L y , L y ] and D Phys = [-L * x , L * x ] × [L * y , L * y ]
(see e.g. Fig. 14 (Left)). Standard PMLs methods [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] require a stretching of the real spatial coordinates following the change of variables

ν = ν + e iθ ν L * ν σ ν (s)ds, ( 58
)
with ν = x, y, and where the absorbing functions are defined by

σ ν (ν) = σ ν (|ν| -L ν ), L * ν |ν| < L ν , 0, |ν| < L * ν .
Finally, we set

S ν (ν) := 1 + e iθν σ(ν),
and define the operators along the ν-direction

∂ γ ν → 1 S γ (ν) ∂ γ ν = 1 1 + e iθ σ(ν) γ ∂ γ ν , ( 59
)
where γ is a given derivation order. Following the 1D case, this leads to various choices of FPMLs for the 2D case by extension along each direction, by adapting γ.

Example 8. Rather than developing the general case (that will be studied in a forthcoming paper dedicated to multidimensional FPDEs), we consider the following FPDE

   i∂ t u(t, x, y) + ∂ 3/2 x u(t, x, y) + 1 5 ∂ 3/2 y u(t, x, y) = 0, t > 0, (x, y) ∈ R 2 , u(t = 0, x, y) = u 0 (x, y), (x, y) ∈ R 2 . ( 60
)
We propose to compute the solution to this system in the bounded domain [0; T ] × D, with D = [-2.5; 2.5] 2 and T = 1. The initial data is chosen with a support close to the boundary (since we are mainly interested in the absorbing features) (see Fig. 14) u 0 (x, y) = e -5((x+5/3) 2 +y 2 )-2ix+iy . Following Approach 1, the modified equation with FPML and periodic boundary conditions reads

       i∂ t u + 1 S 3/4 x ∂ 3/4 x 1 S 3/4 x ∂ 3/4 x u + 1 5 1 S 3/4 y ∂ 3/4 y 1 S 3/4 y ∂ 3/4 y u = 0, in ∈ [0; T ] × D,
+ periodic BC on [0; T ] × ∂D, u(0, x, y) = u 0 (x, y), (x, y) ∈ D.

(61)

We apply a directional splitting with semi-implicit discretization for solving the IBVP. The real space grid has N x × N y = 121 2 points and the time step is ∆t = 5 × 10 -2 . In Fig. 15, we compare snapshots of the amplitude of the solution at times t = 0.1, 0.4, 1 (and their logarithm for the last line at t = 1) for i) a solution of reference (computed on a larger domain), ii) a solution with periodic boundary conditions (without FPML) and iii) a solution with FPML by considering the Type VI profile σ 0 /ν 2 -σ 0 /δ 2 ν , with σ 0 = 10 -3 , θ = π/64, δ ν = 5 × 10 -2 L ν (for L ν = 2.5). As we can observe, the FPML solution reproduces correctly the reference solution. When using the logscale representation, we can see that small residuals of the wave field pass through the bottom boundary and come back at the top interface since we are using periodic boundary conditions. In comparison, the full wave is transmitted from the bottom to the top boundary because of the periodic boundary condition.

The absorbing layer used in the above computations is based on (59) leading to (61) (i.e. using one of the previous 1D strategies). For completeness of the study, we also consider the solution based on the absorption operator built on a direct approach, i.e. replacing the operator in (60) by where α(= 3/2) is the complete derivation order for each direction ν = x, y. To this aim, we report in Fig. 16, the logarithm of the amplitude of the solution at t = 1, by using (62) (Left) and (59) (Right), which confirms that (59) should be rather selected, as discussed in Section 2.

∂ α ν → 1 S α (ν) ∂ α ν = 1 1 + e iθ σ(ν) α ∂ α ν , (62)
u(t, x, y) = 0, t > 0, (x, y) ∈ R 2 , u(0, x, y) = u 0 (x, y), (x, y) ∈ R 2 . ( 63
)
The initial data is chosen as the L 2 -normalized gaussian The discretization scheme uses a directional time splitting with semi-implicit discretization of the IBVP. The spatial grid involves N x × N y = 101 2 points while ∆t = 2.5 × 10 -2 . We report in Fig. 17 a snapshot of the amplitude of the the wave field at t = 0.7 with periodic boundary conditions (without any FPML), the reference solution computed on a large domain, and the FPML solution. Here, we are using the Approach 1 (see Eq. ( 65)) for the Type II profile, setting σ 0 = 5 × 10 -4 , θ = π/4, and δ ν = 10 -1 L ν . A very good absorption is observed with the proposed methodology. Let us remark that some tuning is still necessary to properly choose the FPML parameters, which is problem-dependent as usual. Example 9bis. We now slightly modify system (63) as

   i∂ t u(t, x, y) + 9/10 x u(t, x, y) + 1 10 3/4 y u(t, x, y) = 0, t > 0, (x, y) ∈ R 2 , u(t = 0, x, y) = u 0 (x, y), (x, y) ∈ R 2 , (66)
where the initial data is given by (64). The FPML system related to Approach 1 is

i∂ t u + 1

Focus on FPMLs for FPDE involving the 2D fractional laplacian

As seen in subsection 2.3, accurate FPMLs can be derived for 1D FPDEs involving fractional laplacians. Since FPMLs can be obtained for RL fractional operators, this is also relatively expected for the 1D fractional laplacian since we have some relations between both operators, like e.g. considering the 1D formula, -∞ < x < +∞,

-(-∆) α u(x) = ∂ 2α ∂|x| 2α u(x). (68)
For the 2D case, this is much less clear how to proceed and a relation as the above one does not exist (see e.g. [START_REF] Kwasnicki | Ten equivalent definitions of the fractional laplace operator[END_REF][START_REF] Lischke | What is the fractional laplacian?[END_REF] ). We discuss the question of the extension of FPMLs to 2D FPDEs and prospect a first possible approach. Nevertheless, this example also shows that further developments are still needed.

Let us consider the time-dependent two-dimensional FPDE

   i∂ t u(t, x, y) + α∈R v [α] (x, y)(-) α u(t, x, y) = 0, for (t, x, y) ∈ R * × R 2 , u(t = 0, x, y) = u 0 (x, y), with (x, y) ∈ R 2 . (69)
The Fourier spectral definition of the fractional laplacian is given by

(-) α u = F -1 ((|ξ x | 2 + |ξ y | 2 ) α F(u)(ξ x , ξ y )),
where ξ x (respectively ξ y ) is the Fourier dual variable in direction x (respectively y), F(u) is the two-dimensional Fourier transform of u and F -1 is the associated inverse Fourier transform. For α = 1, let us recall that the modified PML laplacian operator has a total symbol given by

σ(-PML ) = -σ( 1 S x ∂ x ( 1 S x ∂ x ) + 1 S y ∂ y ( 1 S y ∂ y )) = -σ( 1 S 2 x ∂ 2 x + 1 S x ∂ x ( 1 S x )∂ x + 1 S 2 y ∂ 2 y + 1 S y ∂ y ( 1 S y )∂ y ) = 1 S 2 x |ξ x | 2 + i 1 S x ∂ x ( 1 S x )ξ x + 1 S 2 y |ξ y | 2 + i 1 S y ∂ y ( 1 S y )ξ y . (70)
By analogy, we propose the 2D FPML laplacian with variable coefficients based on the pseudodifferential operator definition

(-) α PML u := F -1 (a α (x, y, ξ x , ξ y )F(u)(ξ x , ξ y )),
with the symbol a α (x, y, ξ x , ξ y ) = (-σ( PML )) α .

From a practical point of view, when α is not an integer, it is no longer possible to split the real space and Fourier variables. In the latter case, a direct implementation is very inefficient, and then requires some approximations. A possible approach is however based on the decomposition of the fractional laplacian : α = α-1 . As a consequence, the following approximation is proposed

α PML ≈ α PML,a := ( PML,cst ) α-1 PML , (71)
where

σ( PML,cst ) = -( 1 S 2 x,cst |ξ x | 2 -i S x,cst S 3 x,cst ξ x + 1 S 2 y,cst |ξ y | 2 -i S y,cst S 3 y,cst ξ y ) (72)
and S x,cst , S y,cst , S x,cst , S y,cst are the respective approximations to S x , S y , S x , S y based on a constant (cst) profile. In practice, we solve on [0; T ] × D the approximation of system (69) by

   i∂ t u(t, x, y) + α∈R v [α] (x, y)(-PML,a ) α u(t, x, y) = 0, for (t, x, y) ∈ [0; T ] × D, u(t = 0, x, y) = u 0 (x, y), for (x, y) ∈ D, (73)
with periodic boundary conditions on [0; T ] × ∂D.

Example 10. To illustrate the approach, we consider the following system (with α = 9/10)

   i∂ t u(t, x, y) + e 9πi/10 20 9/10 u(t, x, y) = 0, for (t, x, y) ∈ R * × R 2 , u(t = 0, x, y) = u 0 (x, y), for (x, y) ∈ R 2 , (74)
with u 0 (x, y) = 10e -15((x-2) 2 +(y-16/3) 2 )+5ix-5iy .

The computational domain is D = [-8, 8] 2 and the final time is T = 350. Following our strategy, the modified FPDE with FPML which is solved consists in replacing 9/10 in (74) by

9/10 PML,a = -1/10 PML,cst PML ,
where PML,cst is defined by expression (72). In this example, we simply take S x,cst = S y,cst = 1 + 10 -2 e iθ and therefore S x,cst = S y,cst = 0. We apply a directional splitting with semiimplicit discretization to solve the corresponding IBVP (73). The real space grid involves N x × N y = 101 2 points and we fix ∆t = 5 × 10 -1 . We report in Fig. 19 the amplitude of the initial data, of the solution with periodic boundary conditions (without FPML), the solution of reference (computed on a larger domain) and the FPML-based solution. For this last solution, we use a profile of Type I, with σ 0 = 1 × 10 1 , θ = π/4, and δ ν = 0.15L ν .

A few results in mathematical and numerical analysis

In this section, we develop some analysis of fractional PDEs, and their approximation with FPML using the pseudospectral method.

Some first well-posedness results for the IVP and IBVP

To analyze the stability of the overall method, let us first focus on the well-posedness of 1D linear equations with constant coefficients, i.e. we consider where {v [α] } α∈R is a finite sequence of real or purely imaginary numbers, R is a finite set of positive rational numbers, and u 0 ∈ L 2 (R). By Fourier transform along x, we get the following expression of the solution, for any time t > 0,

   ∂ t u(t, x) + α∈R v [α] ∂ α x u(t, x) = 0, (t, x) ∈ R * × R, u(0, x) = u 0 (x), x ∈ R, (76)
u(t, x) = F -1 x exp -t α∈R v [α] (iξ) α * x u 0 (x),
where * x denotes the convolution product with respect to x. In general, we have u(t,

•) L 2 (R) = u 0 L 2 (R) since R F(u)(t, ξ) 2 dξ = R exp -t α∈R v [α] (iξ) α 2 F(u 0 )(ξ) 2 dξ .
Let us remark that (iξ) α = cos(απ/2)ξ α + i sin(απ/2)ξ α , for α ∈ R and ξ ∈ R. Now, if α = β/γ ∈ Q * + is such that γ ∈ 2N + 1, then ξ α is real, and if γ ∈ 2N * then ξ α = e iαπ/2 |ξ| α if ξ < 0. We then write R = R 1 ∪ R 2 , where R 1 (resp. R 2 ) is composed by relative prime numbers β/γ such that γ is odd (resp. even). This yields

R F(u)(t, ξ) 2 dξ = R α∈R 1 exp -2v [α] (cos(πα/2) + i sin(πα/2))ξ α t × α∈R 2 exp -v [α] (h α (ξ) + ik α (ξ))t 2 F(u 0 )(ξ) 2 dξ,
where

h α (ξ) = cos(απ/2)ξ α , if ξ > 0, cos(απ)|ξ| α , if ξ < 0 , k α (ξ) = sin(απ/2)ξ α , if ξ > 0, sin(απ)|ξ| α , if ξ < 0 .
The well-posedness of the IVP is dependent on R and the sign of {v [α] } αR . Assuming for instance that {v [α] } α∈R are purely complex numbers, then we have

R F(u)(t, ξ) 2 dξ = R α∈R 1 exp -2Im{v [α] } sin(πα/2)ξ α t × α∈R 2 exp -2Im{v [α] }k α (ξ)t F(u 0 )(ξ) 2 dξ . Similarly, if {v [α] } α∈R is real, one gets R F(u)(t, ξ) 2 dξ = R α∈R 1 exp -2v [α] cos(πα/2)ξ α t × α∈R 2 exp -2v [α] h α (ξ)t F(u 0 )(ξ) 2 dξ .
In this case, the well-posedness depends on α and {v [α] } α∈R . For instance, for R ⊂ 2N #R (R 1 = ∅) and {v [α] } α some purely complex numbers, we easily check that we have the norm conservation property u(t,

•) L 2 (R) = u 0 L 2 (R) . If R ⊂ 4N #R
and {v [α] } α is real and positive, then one gets the inequality: u(t, •) L 2 (R) u 0 L 2 (R) . Another example is as follows. Let us e.g. take R = R 2 = {1/4, 1/2} and {v [α] } α positive, then we have

R F(u)(t, ξ) 2 dξ = R α∈R 2 exp -2v [α] h α (ξ)t F(u 0 )(ξ) 2 dξ = R + α∈R 2 exp -2v [α] cos(απ/2)ξ α t F(u 0 )(ξ) 2 dξ + R -α∈R 2 exp -2v [α] cos(απ)|ξ| α t F(u 0 )(ξ) 2 dξ u 0 2 .
In fact, it can be shown that u(t, •) L 2 (R) tends towards zero, when t goes to infinity. More generally, the following proposition holds.

Proposition 5.1. Let us assume that u 0 ∈ L 2 (R) and {v [α] } α∈R are positive or purely complex numbers in iR -, and R is composed by irreducible rational numbers β/γ 1/2 such that γ is even. Then the following IVP

∂ t u(t, x) + α∈R v [α] ∂ α x u(t, x) = 0, u(0, •) = u 0 , is well-posed in C ∞ (I; H r (R))
, where

H r (R) := v ∈ L 2 (R) : F -1 [(1 + |ξ| 2 ) r/2 Fv] ∈ L 2 (R)
and

v H r (R) = R |v| 2 dx + R 2 |v(x) -v(y)| 2 |x -y| 1+2r dxdy 1/2
.

The above results can be extended to inhomogeneous equations with a right hand side f ∈ L 2 (R).

A deeper analysis is necessary for initial boundary value problems on D with FPML. We follow the strategy developed in [START_REF] Li | Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[END_REF], setting I = (0; T ). We refer to [START_REF] Di Nezza | Hitchhiker's guide to the fractional Sobolev spaces[END_REF] for the mathematical study of fractional Sobolev spaces W s,p . In the following, we also define

B r (I × D) := C 0 (I; L 2 (D)) ∩ L 2 (I; H r/2 (D)), (77)
equipped with the norm u B r (I×D) := sup

0 t T v(t, •) L 2 (D) + v L 2 (I;H r/2 (D)) 1/2
.

From [START_REF] Li | Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[END_REF], we recall that, for •, • the L 2 (D)-inner product on D, we have

∂ α x u, v = u, x ∂ α v ,
where x ∂ α is defined in subsection 2.1, and

H α 0 (D) = C ∞ 0 (D) H α (D)
. For α ∈ (0, 2) -{1} and u ∈ H α/2 0 (D), we have

∂ α x u = ∂ α/2 x ∂ α/2 x u .
Let us first prove the existence of solutions for the corresponding stationary problem.

Proposition 5.2. We impose a homogeneous Dirichlet boundary condition u = 0 at ∂D. We assume that the order r of the FPDE is between 1 and 2, i.e. 1 < r < 2. Let us consider the FPDE with FPML based on Approach 1

     α∈R v [α] (x) 1 S βα (x) ∂ βα x nα+1 u(x) = f (x), in D, u(x) = 0, on ∂D, (78)
where {v [α] } α∈R is a finite sequence of complex numbers, u 0 ∈ L 2 (D), f ∈ L 2 (D) and R is a finite set of positive real numbers, where S is assumed to be smooth and defined as in Eq. [START_REF] Lischke | What is the fractional laplacian?[END_REF], with σ 0 or δ x small enough. Then, the BVP (78) is well-posed in H r/2 0 (D). In addition, there exists c > 0 such that the solution satisfies the inequality

u H r/2 0 (D) c f H r/2 0 (D) .
Proof. As n α = 1, we directly have

     α∈R v [α] 1 S βα (x) ∂ βα x 1 S βα (x) ∂ βα x u(x) = f (x), in D, u(x) = 0, on ∂D, (79)
for f ∈ L 2 (D). By taking the L 2 (D)-inner product on D with an element v ∈ H r/2 0 (D), the corresponding variational formulation reads

α∈R v [α] 1 S βα (x) ∂ βα x 1 S βα (x) ∂ βα x u, v = f, v .
By integration by parts, one gets

α∈R v [α] 1 S βα ∂ βα/2 x u, 1 S βα x ∂ βα/2 v + 1 S βα ∂ βα/2 x u, x ∂ βα/2 1 S βα v = f, v .
We set

A(u, v) := α∈R v [α] 1 S βα ∂ βα/2 x u, 1 S βα x ∂ βα/2 v + 1 S βα ∂ βα/2 x u, x ∂ βα/2 1 S βα v . (80)
Let us introduce T [βα] := 1/S βα . Denoting the sup-norm on D by • ∞ , there exists C > 0 such that

1 S βα ∂ βα/2 x u, 1 S βα x ∂ βα/2 v T [βα] 2 ∞ ∂ βα/2 x u L 2 (D) x ∂ βα/2 v L 2 (D) c T [βα] 2 ∞ u H βα/2 (D) v H βα/2 (D) .
The last inequality is a direct consequence of Lemma 2.6 in [START_REF] Li | Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[END_REF], which states that ∂ βα/2 x u L 2 (D) and x ∂ βα/2 v L 2 (D) are equivalent semi-norms on H βα/2 0 (D). Next, by the Cauchy-Schwarz inequality, one gets

1 S βα ∂ βα/2 x u, x ∂ βα/2 1 S βα v x ∂ βα T [βα] ∞ ∂ βα/2 x u L 2 (D) v L 2 (D) x ∂ βα T [βα] ∞ u H βα/2 (D) v H βα/2 (D) .
Therefore, there exists C = C(D; β α ) > 0 such that

|A(u, v)| C α∈R |v [α] | T [βα] 2 ∞ + x ∂ βα T [βα] ∞ u H βα/2 (D) v H βα/2 (D) .
Moreover from Proposition 2.1 in [START_REF] Di Nezza | Hitchhiker's guide to the fractional Sobolev spaces[END_REF], we have

u H βα/2 (D) C u H r/2 (D) ,
which proves, on one hand, the continuity of A in H

(u, v) in H r/2 0 (D) × H r/2 0 (D) |A(u, v)| C u H r (D) v H r (D) .
On the other hand, the coercivity of A involving A(u, u) which writes

A(u, u) := α∈R v [α] 1 S βα ∂ βα/2 x u, 1 S βα x ∂ βα/2 u + 1 S βα ∂ βα/2 x u, x ∂ βα/2 1 S βα u , (81)
is a consequence of Lemma 2.6 from [START_REF] Li | Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[END_REF] and Poincaré's inequality which occurs in H βα/2 0 (D). It is indeed proven in [START_REF] Xu | A note of fractional order Poincaré inequalities[END_REF] that

∂ βα/2 x u 2 L 2 (D) 1 2 (2Vol(D)) -2βα u 2 L 2 (D) .
This implies that, for any α ∈ R, there exists C > 0 such that

1 S βα ∂ βα/2 x u, 1 S βα x ∂ βα/2 u C u 2 H r/2 (D) . (82)
Next, we notice that S βα is constant outside the layer, leading to

Re 1 S βα ∂ βα/2 x u, x ∂ βα/2 1 S βα u = Re [-L,-L * ]∪[L * ,L] 1 S βα (x) x ∂ βα/2 1 S βα (x) u(x)∂ βα/2 x u(x)dx .
Furthermore, the term x ∂ βα/2 1/S βα is linear in σ 0 . Now, we easily deduce that there exists

c(σ 0 , δ x , L) > 0 such that Re 1 S βα ∂ βα/2 x u, x ∂ βα/2 1 S βα u -c(σ 0 , δ x , L) u 2 H βα/2 (D) . (83)
Finally, we conclude that, by taking δ x (= L -L * ) or σ 0 small enough, (83) can be controlled by (82). As a consequence, there exists

C = C(σ 0 , δ x , L) > 0 such that ReA(u, u) C u 2 H r/2 (D) .
Now, the continuity of L : v ∈ H r/2 0 (D) → f, v is obvious from Cauchy-Schwarz's inequality and since f ∈ L 2 (D). From the Lax-Milgram theorem [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], we finally prove the existence of solution in H r/2 (D) to (79), ending hence the proof.

Before considering the full time-dependent problem, we recall that H r/2 (D) is a Hilbert space and that the embedding of H r/2 (D) to L 2 (D) is compact. Assuming now that the sequence {v [α] } α has a negative real part and taking A defined in (81) as a symmetric, continuous and elliptic operator, there exists an hilbertian basis {w m } m in H r/2 (D) such that one gets the spectral decomposition

A(w m , v) = λ m w m , v , ∀v ∈ H r/2 0 (D)
, where λ m is an increasing sequence of bounded eigenvalues. This allows us to prove the following Proposition.

Proposition 5.3. Let us assume that the order r of the FPDE is between 1 and 2, i.e. 1 < r < 2. We consider the modified FPDE with homogeneous Dirichlet boundary conditions and FPML related to Approach 1, that is

         ∂ t u(t, x) + α∈R v [α] 1 S βα (x) ∂ βα x nα+1 u(t, x) = f (t, x), (t, x) ∈ R * × D, u(t = 0, x) = u 0 (x), x ∈ D, u(t, x) = 0, (t, x) ∈ R * × ∂D, (84)
where {v [α] } α∈R is a finite sequence of negative real valued numbers1 , u 0 ∈ L 2 (D), f ∈ L 2 (I; L 2 (D)) and R is a finite set of positive rational numbers. The function S βα is assumed to be smooth and defined as in Eq. (34), with σ 0 or δ x small enough. Then, this problem is well-posed in B r (I × D) defined by (77).

Sketch of the proof. We mainly follow the steps presented in [START_REF] Raviart | Introduction à l'Analyse Numérique des Equations Dérivées Partielles[END_REF]. As n α = 1, we first have

∂ t u(t, x) + α∈R v [α] 1 S βα (x) ∂ βα x 1 S βα (x) ∂ βα x u(t, x) = f, u(0, •) = u 0 . For v ∈ H r/2 0 (D), we consider the full variational formulation d dt u(t, •), v + α∈R v [α] 1 S βα (x) ∂ βα x 1 S βα (x) ∂ βα x u, v = f, v .
By integration by parts (see Lemma 2.4 [START_REF] Li | Existence and uniqueness of the weak solution of the space-time fractional diffusion equation and a spectral method approximation[END_REF]), one gets

d dt u(t, •), v + α∈R v [α] 1 S βα ∂ βα/2 x u, 1 S βα x ∂ βα/2 v + 1 S βα ∂ βα/2 x u, x ∂ βα/2 1 S βα v = f, v .
We then seach for u ∈ L 2 (I; H r/2 0 D) ∩ C 0 (I; L 2 (D)). From the existence of a hilbertian basis, we can introduce a sequence {u m } m : t → u m (t, •) which is a Cauchy sequence in L 2 (I; H The generalization to order r-FPDEs, with r > 0 and r ∈ 2N * , is not trivial. We search for u ∈ B r (I × D) such that, for any v ∈ H r/2

0 (D), d dt u(t, •), v + α∈R v [α] 1 S βα (x) ∂ βα x nα+1 2 u, 1 S βα (x) ∂ βα x nα+1 2 v = f, v .
By using similar arguments, we conjecture that the problem is still well-posed but have no rigorous mathematical proof yet.

Some partial stability results of the discretization schemes

In this subsection, we state a stability result. A semi-explicit splitting scheme approximating [START_REF] Xu | A note of fractional order Poincaré inequalities[END_REF] reads

               u n+1/2 h = u n h + ∆t 2 f x h , u n+1/2 h , u n * h = u n+1/2 h -∆t α∈R v [α] h 1 S βα h [[∂ βα x ]] nα+1 u n+1/2 h , u n+1 h = u n * h + ∆t 2 f x h , u n+1 h , (85)
while an implicit one is given by system (51). An alternative method is the following semiimplicit Crank-Nicolson scheme

               u n+1/2 h = u n h + ∆t 2 f x h , u n+1/2 h , I + ∆t 2 α∈R v [α] h 1 S βα h [[∂ βα x ]] nα+1 u n * h = I - ∆t 2 α∈R v [α] h 1 S βα h [[∂ βα x ]] nα+1 u n+1/2 h , u n+1 h = u n * h + ∆t 2 f x h , u n+1 h . (86)
The L 2 -stability of (85), given by

u n+1/2 h 2 I + ∆t α∈R v [α] h 1 S βα h [[∂ βα x ]] nα+1 u n * h 2 ,
requires to state the estimates

ρ I -∆t α∈R v [α] h S [qα] h [[∂ qα x ]] 1 S [1] h [[∂ x ]] nα 1 ,
while for (51) one needs to prove

ρ I + ∆t α∈R v [α] h S [qα] h [[∂ qα x ]] 1 S [1] h [[∂ x ]] nα -1
1 .

The semi-implicit Crank-Nicolson scheme (86) is stable if

ρ I - ∆t 2 α∈R v [α] h 1 S βα h [[∂ βα x ]] nα+1 I + ∆t 2 α∈R v [α] h 1 S βα h [[∂ βα x ]] nα+1 -1 1 .
Computing these various radii is clearly far from being trivial. In the case of the explicit scheme (85), we conjecture that standard explicit real space methods for solving linear FPDE of the form u t + γ∂ α x u = 0, with α > 0 and γ ∈ R * , are stable under a CFL condition of the form [START_REF] Baeumer | Boundary conditions for fractional diffusion[END_REF][START_REF] Guo | Fractional Partial Differential Equations and Their Numerical Solutions[END_REF][START_REF] Karniadakis | Special Issue on Fractional PDEs: Theory, Numerics, and Applications[END_REF][START_REF] Lischke | What is the fractional laplacian?[END_REF] ∆t ∆x α C(γ), for some C(γ) > 0. For α ∈ N * and on unbounded grids, this is a consequence of D α x := max ξ∈[-π/h,π/h] |iξ| α = (π/h) α , D α

x being an order α-derivative operator. For the periodic case, and an even number N x of points, we have: if α is an even integer, D α Nx = π/h -1 α , while when α is an odd integer, D α Nx = π/h α .

To analyze the method, we follow the strategy presented in [START_REF] Kreiss | Stability of the Fourier method[END_REF] that we apply to FPDEs without FPML.

where {v [α] } α∈R are real (resp. purely imaginary) smooth functions and R is a finite set of positive real numbers such that α ∈ R, cos(α/2)v [α] (x) 0 (resp. sin(α/2)v [α] (x) 0) for any α ∈ R and x ∈ [-L, L], such that (87) is well-posed. Then, the implicit scheme

u n+1 h + ∆t α∈R v [α] h [[∂ α x ]]u n+1 h = u n h is unconditionally stable.
Let us remark that the assumption cos(α/2)v [α] (x) 0, for any α ∈ R and x ∈ [-L, L], is not optimal and could be improved by additional technical computations. Next, let us recall that, for any function w w(x) = h u h . As in [START_REF] Kreiss | Stability of the Fourier method[END_REF], T h is a skew-hermitian matrix, so that

T [α] h h := sup u h h =1 T [α] h u h h = (πN x ) α ,
and the eigenvalues of T h u n+1 h .

Let us define Λ

[α]

h = diag λ [α] 0 , • • • , λ [α]
Nx-1 . We set

V [α] h = diag v [α] (x 0 ), • • • , v [α] (x Nx-1
) ∈ M Nx (C), 0 h ∈ C Nx , I h as the identity matrix in C Nx×Nx , and

du h dt + α∈R V [α]
h T

[α]

h u h = 0 h which is approximated by

u n+1 h -u n h ∆t + α∈R V [α] h T [α] h u n+1 h = 0 h ,
or equivalently

I h -∆t α∈R V [α] h T [α] h u n+1 h = I h u n h .
We introduce the norm • h defined by

u 2 h = u, u h = h Nx-1 j=0
u(x j )u(x j ).

Let us recall that e 2inx , e 2imx h = 1 if m = n and 0 if |m -n| N x . We denote by Π h the matrix of eigenvectors to T

[α] h = Π h Λ [α] h Π -1
h , and by w h := Π -1 h u h , so that

I h -∆t α∈R V [α] h Λ [α] h w n+1 h = I h w n h . As V [α] h Λ [α]
h is diagonal, we get, from λ w n h (x j )w n+1 h (x j )

h|w n+1 h (x j )| 2 1/2 h |w n h (x j )| 2 1 -λ [α] j v [α] (x j )∆t 2 1/2
. Now, if we assume that v [α] is real then we obtain

1 -λ [α] j v [α] (x j )∆t| = (1 -2(2πj) α v [α] (x j )∆t cos(α/2)Γ(2)/Γ(2 -α) + (2πj) 2α (v [α] (x j )) 2 ∆t 2 Γ 2 (2)/Γ 2 (2 -α)
Since v [α] (x j ) cos(α/2) 0, this yields, for any n 0,

w n+1 h w n h .
The latter occurs similarly if v [α] is purely complex and v [α] (x j ) sin(α/2) 0. Finally, one gets

u n+1 h u n h ,
concluding then the proof.

A similar result can be stated for the semi-implicit scheme (86).

where κ ∈ R, {v [α] } α∈R are real (resp. purely imaginary) smooth functions and R is a finite set of positive rational numbers such that for any α ∈ R, cos(α/2)v [α] 0 (resp. sin(α/2)v [α] 0). Then, the scheme

           u n+1/2 h = u n h + iκ ∆t 2 u n h 2 u n+1/2 h , u n+1 h = u n+1/2 h -∆t α∈R v [α] h [[∂ α x ]]u n+1 h , u n+1 h = u n * h + iκ ∆t 2 u n * h 2 u n+1 h , (89)
is L 2 unconditionally stable.

Proof. We trivially have that

u n+1/2 h u n h , u n+1 h u n * h ,
which combined with Proposition 5.4 leads to the unconditional stability of the scheme (89).

Corollary 5.2. We consider now the FPDE

∂ t u(t, x) + α∈R v [α] (x) S α (x) ∂ α x u(t, x) = 0, ( 90
)
where {v [α] } α∈R are real (resp. purely imaginary) smooth functions and R is a finite set of positive rational numbers such that α ∈ R, cos(α/2)v [α] 0 (resp. sin(α/2)v [α] 0) for any α ∈ R and r = max R < 1. The implicit scheme

u n+1 h + ∆t α∈R v [α] h S α h [[∂ α x ]]u n+1 h = u n h
is unconditionally stable for the L 2 -norm.

Conclusion

The aim of this paper is to introduce some possible extensions of the PML method to FPDEs. We first start with one-dimensional problems by building fractional PMLs (FPMLs) and then various approximations are derived. Next, we propose an implementation of these FPMLs by using pseudospectral techniques in space variables and time splitting schemes. Some numerical examples show that the resulting scheme with FPMLs is efficient and accurate. Extensions to FPDEs involving fractional laplacians and to higher-dimensional problems are developed. Finally, some mathematical preliminary results are stated about the well-posedness of the truncated problems. In addition, a numerical analysis of the schemes is proposed.

Figure 1 :

 1 Figure 1: Example 1. Comparison of the amplitude |u| (logscale) for the Approaches 0, 1, 2a, 2b, 3, with the FPML of Type I, pure Periodic Boundary Conditions (i.e. without any PML), and a reference solution computed on a large domain.

Figure 2 :

 2 Figure 2: Example 2. Comparison of the amplitude |u| (logscale) for the Approaches 0, 1, 2a, 2b, 3, with the FPML of Type I, pure Periodic Boundary Conditions (i.e. without any PML), and reference solution computed on a large domain.

Figure 3 :

 3 Figure 3: Example 3. Solutions {(x, t, log |u(x, t)|), (x, t) ∈ D × [0; T ]} with periodic boundary conditions (Left), FPML of Type I (Center) and reference solution (Right).

Figure 4 :

 4 Figure 4: Example 4. Solution {(x, t, log |u(x, t)|), (x, t) ∈ D × [0; T ]} with periodic boundary conditions (Left), FPML of Type I (Center) and reference solution (Right).

Figure 5 :

 5 Figure 5: Example 4. Comparison the FPML with Type I, II and VI absorption profiles.

Figure 6 :

 6 Figure 6: Example 4. Comparison of the FPML of Type I with different values of σ 0 (= 10 -1 , 1, 10, 10 2 , 10 3 ).

  (Left) the absorption effects of the FPML (Approach 2a) for three different profiles: Type I-, Type II (with σ 0 = 10) and Type VI (with σ 0 = 10 -1 ), setting δ x = 0.05L and θ = π/8. In addition, we also report on Fig. 9 (Right) the L 2 (D)-norm of the solution (in logscale) vs. the time variable t, i.e. {(t, log u(t, •) L 2 (D) ) : t ∈ [0; T ]}.

Figure 7 :

 7 Figure 7: Example 5. Solution {(x, t, log |u(x, t)|), (x, t) ∈ D × [0; T ]} with periodic boundary conditions (no FPMLs) (Left), with FPML of Type I (Center) and reference solution (Right).

Example 6 .

 6 Let us now consider the following nonlinear cubic Schrödinger equation

Figure 8 :

 8 Figure 8: Example 5. Solution {(x, t, log |u(x, t)|), (x, t) ∈ D × [T /2, T ]} with periodic boundary conditions (Left), FPML with absorption function of Type I (Center), and reference solution (Right).

Figure 9 :

 9 Figure 9: Example 5. (Left) Amplitude |u| of the solution (in logscale) at final time T = 7.5, on the domain D, for the FPMLs (Approach 2a) with absorption function of Type I, II (σ 0 = 10) and VI (σ 0 = 10 -1 ), for θ = π/8. (Right) Time evolution of the L 2 (D)-norm of the solution for the FPMLs with absorption function of Types I, II, VI, with periodic boundary conditions and for the reference solution.

Figure 10 :

 10 Figure 10: Example 6. Solution {(x, t, log |u(x, t)|), (x, t) ∈ D × [0; T ]} with periodic boundary conditions (Left), FPML with absorption function of Type I (Center) and reference solution (Right).

Figure 11 :

 11 Figure 11: Example 6. log |u(t, x)| at time T = 1.25 with periodic boundary conditions, FPML with Type I absorption function for σ 0 = 5 × 10 1 , δ x = 5 × 10 -2 L and θ = π/8.

Figure 12 :

 12 Figure 12: Example 7. Solution {(x, t, log |u(x, t)|), (x, t) ∈ D × [0; T ]} with periodic boundary conditions (Left), FPML with absorption function of Type I (Center) and reference solution (Right).

Figure 13 :

 13 Figure 13: Example 7. (Left) Amplitude |u| (logscale) of the solution at final time T = 10 and on D = [-8, 8]. The FPMLs (Approach 2a) are based on the absorption functions of Type I, II (σ 0 = 10) and VI σ 0 = 10 -1 , with θ = π/8. (Right) Amplitude |u| (logscale) of the solution at final time T = 10 on D = [-10, 10]. We use the FPML with absorption profiles of Types I, II, θ = π/8, σ 0 = 10 2 and for the Approaches 1 and 2.

Figure 14 :

 14 Figure 14: Example 8. (Left) Domain with FPML. (Right) Amplitude of the initial data.

Figure 15 :

 15 Figure 15: Example 8. Three first lines (for the times t = 0.1, 0.4, 1): amplitude of the solution: (Left) reference, (Center) FPML with Type VI profile, (Right) periodic BC. Fourth line: logarithm of the amplitude of the solution at time t = 1: (Left) reference, (Center) FPML with Type VI profile, (Right) periodic BC.

Figure 16 :

 16 Figure 16: Example 8. Logarithm of the amplitude of the solution at time t = 1: (Left) by using the direct formulation (62), (Right) by using (59).

u 0 1 S

 01 (x, y) = N -1 e -5((x+10/3) 2 +y 2 )-10ix-2iy . (64)We consider the bounded spatial domain of computation D = [-5, 5] 2 while the time interval is [0; T ], with T = 1. The truncated FPDE with FPML is then based on Approach 1 following i∂ t u + T ]×D, adding periodic boundary conditions at ∂D and considering the initial data (64).

Figure 17 :

 17 Figure 17: Example 9. (Top-left) amplitude of the initial data, (Top-right) periodic solution, (Bottom-left) reference solution and (Bottom-right) of the FPML solution with Type II profile.

  T ] × D, with D = [-5, 5] 2 and for T = 1. The scheme as well as the discretization parameters are the same as for Example 9. We report in Fig.18the amplitude of the initial data, and the solution (at time t = 0.63) with periodic boundary conditions (no FPML), a reference solution, and the FPML solution based on the Type II profile, with σ 0 = 4 × 10 -4 , θ = π/4, and δ ν = 10 -1 L ν . We again clearly observe an effective accuracy of the FPML implemented within the pseudospectral scheme.

Figure 18 :

 18 Figure 18: Example 9bis. (Top-left) amplitude of the initial data, (Top-right) periodic solution, (Bottomleft) reference solution and (Bottom-right) of the FPML solution with Type II profile.

Figure 19 :

 19 Figure 19: Example 10. (Top-left) amplitude of the initial data, (Top-right) periodic solution, (Bottom-left) reference solution and (Bottom-right) of the FPML solution with Type I profile.

0

  (D) since there exists C > 0 such that for all

r/ 2 0 2 0 2 0 2 0

 2222 (D)) ∩ C 0 (I; L 2 (D)), both being complete spaces. The continuous embedding of L 2 (I; H r/(D)) and C 0 (I; L 2 (D))) in L 2 (I; L 2 (D)) justifies that the limit of {u m } m is identical in L 2 (I; H r/(D)) and C 0 (I; L 2 (D)). Finally, since C ∞ 0 (D) is dense in H r/(D), we can conclude.

Proposition 5 . 4 .∂

 54 We consider the FPDE t u(t, x)+ α∈R v [α] (x)∂ α x u(t, x) = 0, in [0; T ] × [-L, L], u(t, -L) = u(t, L) t > 0, u(0, x) = u 0 (x), in [-L, L],

Nx/ 2 k=-Nx/ 2 a 2 (

 222 (k)e 2πix , then, for {x j } 0 j Nx-1 , we have the representation 2iπk) α a(k)e 2πikx j .Proof. Let uh = u h (x 0 ), • • • , u(x Nx-1 )T and y[α] h = ∂ α x u(x 0 ), • • • , ∂ α x u(x Nx-1 )

h

  are given by λ[α] k = (2πik) α Γ(2)/Γ(2 -α) with the eigenvectors 1, e 2πikhx , • • • , e 2πikNxhx with k = 0, • • • N x .As a consequence, we can rewrite α∈R v

  ) α = (2πj) α e iα/2 Γ(2)/Γ(2 -α), the inequality h j v [α] (x j )∆t

Corollary 5 . 1 .

 51 Let us consider the nonlinear FPDE∂ t u(t, x) + α∈R v [α] (x)∂ α x u(t, x) = iκ|u| 2 u(t, x),

The result could be extended to complex numbers {v[α] } α∈R with negative real part.

Acknowledgements. X. Antoine was supported by the French National Research Agency project NABUCO, grant ANR-17-CE40-0025. E. Lorin was supported by the Discovery Grant program of the Natural Sciences and Engineering Research Council of Canada (NSERC).