Christophe Chesneau
email: christophe.chesneau@unicaen.fr

Hassan S Bakouch
email: hnbakouch@yahoo.com

Nauman Muhammad

Khan

A weighted transmuted exponential distributions with environmental applications

Keywords: Conditional moments, Estimation, Statistical distributions, Generating function, Return period. 2000 MSC: 60E05, 62E15

In this paper, we introduce a new three-parameter distribution. It is based on the combination of a re-parametrization of a general family of distributions (known as the EGNB2 distribution) and the so-called quadratic rank transmutation map defined with the exponential distribution as baseline. We explore some mathematical properties of this distribution including the hazard rate function, moments, the moment generating function, the quantile function, various entropy measures and (reversed) residual life functions. A statistical study investigates estimation of the parameters using the method of maximum likelihood. The distribution along with other existing distributions are fitted to two environmental data sets and its superior performance is assessed by using some goodness-of-fit tests. As a result, some environmental measures associated with these data are obtained such as the return level and mean deviation about this level.

Introduction

The precise analysis of a wide variety of data sets is limited by the use of models based on the classical distributions (normal, exponential, logistic. . . ). For instance, the analysis of environmental data sets collecting from observations of complex natural phenomena needs special treatments to reveal all the underlying informations. Over the last decades, numerous solutions have been provided by the statisticians, including the elaboration of several methods which aim to increase the flexibility of the former classical distributions. Among these methods, a popular one is to construct generator of distributions by compounding continuous distributions with well-known discrete distributions. This compounding is always motivated by practical problems as those involving cdf of minimum or maximum of several independent and identically random variables. An exhaustive survey on the construction of such generators, with the presentation of new ones, can be found in [START_REF] Tahir | Compounding of distributions: a survey and new generalized classes[END_REF], and the references therein. Among the long list, let us briefly present the EGNB2 distribution introduced by [20, Remark 2 (ii)]. Using a cumulative distribution function (cdf) G(x), the general form of the associated cdf is given by

F EGN B2 (x) = [1 + η υG(x) α ] -1 η -1 (1 + η υ) -1 η -1 . (1)
The EGNB2 distribution can be viewed as an extension of the G-negative binomial families introduced by [START_REF] Louzada | The exponential negative-binomial distribution: A continuous bridge between under and over dispersion on a lifetime modelling structure[END_REF] and [START_REF] Percontini | The G-negative binomial family: General properties and applications[END_REF]. It enjoys remarkable theoretical and practical properties.

In this study, we consider a particular case of this EGNB2 distribution consisting in a re-parametrization for the parameters α, η and υ (1) as described below. Let γ > 0, η = -γ γ+1 , υ = -(γ + 1) and α = 1. That yields a cdf of the (simple) form:

F (x) = [1 + γ G(x)] 1 γ +1 -1 (1 + γ) 1 γ +1 -1 . (2)
Let us now explain the importance of this re-parametrization of [START_REF] Ahammed | Variability of annual daily maximum rainfall of Dhaka, Bangladesh[END_REF], with some statistical features. One can observe that F (x) as the following integral form: F (x) = G(x) 0 p(t) dt, where p(t) denotes the pdf: p(t) = γ+1 (1+γ)

1 γ +1 -1
(1 + γ t) 1 γ . So it reveals to be a new particular case of the T-X family cdf introduced by [START_REF] Alzaatreh | A new method for generating families of continuous distributions[END_REF]. Another remark is that, when G(x) → 0, we have F (x) ∼ γ+1 (1+γ)

1 γ +1 -1
G(x) and when γ → 0, we have F (x) ∼ e G(x) -1 e-1 . This transformation of cdf corresponds to the one proposed in [START_REF] Kumar | A Method of Proposing New Distribution and its Application to Bladder Cancer Patients Data[END_REF]. All the resulting distributions have demonstrated nice properties in terms of analysis of real life data sets. Furthermore, let us observe that the probability density function (pdf) associated to [START_REF] Akinsete | The beta-Pareto distribution[END_REF] is given by

f (x) = (γ + 1) [1 + γ G(x)] 1 γ g(x)
(1 + γ)

1 γ +1 -1
.

Note that we can also express it as a weighted pdf: f (x) = c w(x) g(x), where w(x) = [1 + γ G(x)]

1 γ
is a weight function and c = γ+1 (1+γ)

1 γ +1 -1
is a normalizing constant. It thus belongs to the family of weighted distributions. Further details on such family of distributions can be found in [START_REF] Saghir | Weighted Distributions: A Brief Review, Perspective and Characterizations[END_REF].

On the other side, a remarkable lifetime distribution was introduced by [START_REF] Aryal | transmuted Weibull distribution: a generalization of the Weibull probability distribution[END_REF].

It is characterized by the cdf G(x) = (θ + 1)H(x) -θ[H(x)] 2 , θ ∈ [-1, 1]
, where H(x) denotes the cdf of an exponential distribution with parameter λ. Among others, it demonstrates a clear superiority in terms of fit in comparison to the former exponential distribution.

In this paper, we introduce a new three-parameter distribution which combines the features of the distribution characterized by [START_REF] Akinsete | The beta-Pareto distribution[END_REF] and the transmuted exponential distribution. We thus obtain a very flexible distribution, which opens new perspectives in terms of the construction of statistical models for data analysis. The theoretical and practical aspects are explored in an exhaustive way. The theoretical ones includes expansions of the cdf, pdf, hrf, quantile function, moments, moment generating function, various entropy measures, residual life functions, conditional moments, mean deviations and reversed residual life function. We investigate the estimation of its parameters via the maximum likelihood method. Two real data sets in environmental sciences are analyzes to show its superior performance in terms of fit in comparison to well-known distributions: the Gamma distribution, the Marshal-Olkin exponential distribution [START_REF] Marshall | A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families[END_REF], the Nadarajah-Haghighi exponential distribution [START_REF] Nadarajah | An extension of the exponential distribution[END_REF], the exponentiated exponential distribution [START_REF] Gupta | Generalized Exponential Distributions[END_REF], the transmuted Weibull distribution [START_REF] Aryal | transmuted Weibull distribution: a generalization of the Weibull probability distribution[END_REF], the transmuted generalized exponential distribution [START_REF] Khan | Transmuted generalized exponential distribution: A generalization of the exponential distribution with applications to survival data[END_REF], the transmuted linear exponential distribution [START_REF] Tiana | Transmuted Linear Exponential Distribution: A New Generalization of the Linear Exponential Distribution[END_REF] and the Kappa distribution [START_REF] Mielke | Another Family of Distributions for Describing and Analyzing Precipitation Data[END_REF]. Then, using our new distribution, estimations of important hydrologic parameters are given.

The rest of this article is organized as follows. In Section 2, we present our main distribution. Some of its mathematical properties are studied in Section 3. Residual life functions are determined in 4. Estimations of the parameters are investigated in Section 5. Applications to two real data sets are provide in Section 6.

A new weighted transmuted exponential distribution

In this section, we precise what is the considered cdf G(x) is given by [START_REF] Akinsete | The beta-Pareto distribution[END_REF]. [START_REF] Shaw | The alchemy of probability distributions: beyond Gram-Charlier expansions, and a skew-kurtotic-normal distribution from a rank transmutation map[END_REF] and [START_REF] Aryal | transmuted Weibull distribution: a generalization of the Weibull probability distribution[END_REF] introduced the quadratic rank transmutation map (QRTM) to propose a new distribution based on the Weibull/exponential one with great flexibility and nice fit for real life data. For these reasons, we use it in our study. We consider the cdf:

G(x) = (θ + 1)H(x) -θ [H(x)] 2 , θ ∈ [-1, 1],
where H(x) is the cdf of the exponential distribution of parameter λ:

G(x) = (θ + 1)(1 -e -λ x ) -θ(1 -e -λ x ) 2 , x, λ > 0.
Putting the above expression into (2), we introduce a new cdf defined by

F (x) = 1 + γ (θ + 1) (1 -e -λ x ) -γ θ (1 -e -λ x ) 2 1 γ +1 -1 (1 + γ) 1 γ +1 -1 = 1 + γ -γ e -λ x (1 -θ + θ e -λ x ) 1 γ +1 -1 (1 + γ) 1 γ +1 -1 , x > 0, λ, γ > 0.
Another useful expression is the following:

F (x) = (1 + γ) 1 γ +1 1 -γ 1+γ e -λ x (1 -θ + θe -λ x ) 1 γ +1 -1 (1 + γ) 1 γ +1 -1 . (3)
We will refer to the distribution given by (3) as the new weighted transmuted exponential and denote it by NWTE(λ, γ, θ) with the considered parameters. The associated pdf is given by

f (x) = λ (1 + γ) 1 γ +1
(1 + γ)

1 γ +1 -1 1 - γ 1 + γ e -λ x 1 -θ + θ e -λ x 1 γ e -λ x 1 -θ + 2 θ e -λ x . (4)
The associated hrf is given by

h(x) = λ 1 -γ 1+γ e -λ x (1 -θ + θe -λ x ) 1 γ e -λ x 1 -θ + 2 θ e -λ x 1 -1 -γ 1+γ e -λ x (1 -θ + θ e -λ x ) 1 γ +1 . (5)
Let us now discuss the possible shapes of pdf (4) and hrf (5) as follows.

lim x→0 f (x) = λ(1 + θ) (1 + γ) (1 + γ) 1 γ +1 -1 , lim x→+∞ f (x) = 0.
On the other hand, we have

lim x→0 h(x) = λ(1 + θ) (1 + γ) (1 + γ) 1 γ +1 -1 , lim x→+∞ h(x) = λ.
In order to visualize the wide variety of shapes, some plots of the pdf (4) and hrf (5) are given below.

Structural properties of the NWTE distribution

Expansion for the associated functions

Expansion for the cdf function. First of all, set h(u) = u(1 -θ + θu), u ∈ (0, 1), θ ∈ [-1, 1]. Note that we have h (u) = 1 -θ + 2θu ≥ min(1 -θ, 1 + θ) ≥ 0, so h is increasing. Since h(0) = 0 and h(1) = 1, we have 0 < h(u) < 1 for all u ∈ (0, 1). Since 0 ≤ γ 1+γ < 1 and 0 < e -λ x (1 -θ + θ e -λ x ) = h(e -λ x ) < 1, the generalized binomial expansion, we have where

1 - γ 1 + γ e -λ x (1 -θ + θ e -λ x ) 1 γ +1 = +∞ i=0 1/γ + 1 i - γ 1 + γ i e -λ i x (1 -θ + θ e -λ x ) i = +∞ i=0 i k=0 1/γ + 1 i i k - γ 1 + γ i e -λix (1 -θ) i-k θ k e -λkx = +∞ i=0 i k=0 H i,k e -λ x (i+k) , (6)
H i,k = 1/γ + 1 i i k γθ -γ 1 + γ i θ 1 -θ k .
Therefore we can expand the cdf function as

F (x) = 1 (1 + γ) 1 γ +1 -1 (1 + γ) 1 γ +1 +∞ i=0 i k=0 H i,k e -λ x (i+k) -1 . (7)
Expansion for the pdf function. Similar mathematical arguments used for [START_REF] Gupta | Generalized Exponential Distributions[END_REF] give

1 - γ 1 + γ e -λ x (1 -θ + θe -λ x ) 1 γ = +∞ i=0 i k=0 A i,k e -λ x (i+k) ,
where

A i,k = 1/γ i i k γ θ -γ 1 + γ i θ 1 -θ k . Therefore f (x) = λ (1 + γ) 1 γ +1
(1 + γ)

1 γ +1 -1 +∞ i=0 i k=0 A i,k e -λ x (i+k) e -λ x 1 -θ + 2 θ e -λ x = +∞ i=0 i k=0 B i,k (1 -θ)e -λ x(i+k+1) + 2 θ e -λ x (i+k+2) , (8)
where

B i,k = λ (1 + γ) 1 γ +1 (1 + γ) 1 γ +1 -1 A i,k .
On the survival function. Note that

S(x) = 1 -F (x) = (1 + γ) 1 γ +1 1 -1 -γ 1+γ e -λ x (1 -θ + θe -λ x ) 1 γ +1
(1 + γ)

1 γ +1 -1 . (9)
Using ( 7), we have the following expansion

S(x) = (1 + γ) 1 γ +1
(1 + γ)

1 γ +1 -1 1 - +∞ i=0 i k=0 H i,k e -λ x (i+k) . ( 10
)
Expansion for the hrf function. Using ( 5), ( 8) and ( 10), an expansion of the hrf function is given by

h(x) = λ +∞ i=0 i k=0 A i,k (1 -θ) e -λ x (i+k+1) + 2 θ e -λx(i+k+2) 1 - +∞ i=0 i k=0 H i,k e -λx(i+k) . ( 11
)
Another expansion comes from the geometric series decomposition:

1 1 -1 -γ 1+γ e -λ x (1 -θ + θe -λ x ) 1 γ +1 = +∞ m=0 1 - γ 1 + γ e -λ x (1 -θ + θ e -λ x ) m( 1 γ +1)
.

By ( 5) and similar mathematical arguments used for [START_REF] Gupta | Generalized Exponential Distributions[END_REF] give:

h(x) = λ +∞ m=0 1 - γ 1 + γ e -λ x (1 -θ + θ e -λ x ) m( 1 γ +1)+ 1 γ (1 -θ) e -λ x + 2 θ e -2λx = +∞ m=0 +∞ i=0 i k=0 G i,k,m (1 -θ) e -λ x (i+k+1) + 2 θ e -λ x (i+k+2) ,
where

G i,k,m = λ m 1 γ + 1 + 1 γ i i k γ θ -γ 1 + γ i θ 1 -θ k .

Quantile function

The quantile functions are in widespread use in general statistics to obtain mathematical properties of a distribution and often find representations in terms of lookup tables for key percentiles. For generating data from the NWTE model, let u ∼ U (0 , 1 ). Then, by inverting the cdf (3) and after some algebra, we get the quantile function

Q(u) = 1 λ       -ln            1 -       1 + θ -(1 + θ) 2 -4 θ γ u [1 + γ] 1 γ +1 -1 + 1 γ γ+1 -1 2 θ                        . ( 12
)
The analysis of the variability of the skewness and kurtosis of X can be investigated based on quantile measures. The Bowley skewness is given by

S = Q(3/4) -2Q(1/2) + Q(1/4) Q(3/4) -Q(1/4)
and the Moors' kurtosis by

K = {Q(7/8) -Q(5/8)} + {Q(3/8) -Q(1/8)} Q(6/8) -Q(2/8) ,
where Q(u) is given by ( 12). These measures are less sensitive to outliers and they exist even for distributions without moments. Figure 3 displays plots of S and K as functions of θ and γ, which show their variability in terms of the shape parameters.

Moments and moment generating function

Moments. Using equation ( 8) and the Gamma function Γ(ν) = +∞ 0

x ν-1 e -x dx, the r-th moments about the origin is given by

µ r = E(X r ) = +∞ -∞ x r f (x)dx = +∞ i=0 i k=0 B i,k (1 -θ) +∞ 0 x r e -λ x (i+k+1) dx + 2 θ +∞ 0 x r e -λ x (i+k+2) dx = Γ(r + 1) λ r+1 +∞ i=0 i k=0 B i,k 1 -θ (i + k + 1) r+1 + 2 θ (i + k + 2) r+1 . ( 13
)
The moment generating function. Similarly the moment generating function associated to the NWTE distribution is given by, for t ≤ λ,

M X (t) = E(e tx ) = +∞ -∞ e t x f (x)dx = +∞ i=0 i k=0 B i,k (1 -θ) +∞ 0 e t x e -λ x (i+k+1) dx + 2 θ +∞ 0 e t x e -λ x (i+k+2) dx = +∞ i=0 i k=0 B i,k 1 -θ λ(i + k + 1) -t + 2 θ λ(i + k + 2) -t . (14)

Entropies

An entropy can be considered as a measure of uncertainty of probability distribution of a random variable. Therefore, we obtain three entropies for the NWTE distribution with investigating a numerical study among them. Entropy 1. Let us consider the Shannon entropy [START_REF] Shannon | Mathematical theory of communication[END_REF]:

H(f ) = -E [log[f (X)]] = - +∞ -∞ f (x) log[f (x)]dx. One can observe that H(f ) = -log λ (1 + γ) 1 γ +1 (1 + γ) 1 γ +1 -1 - 1 γ +∞ 0 f (x) log 1 - γ 1 + γ e -λ x (1 -θ + θ e -λ x ) dx + λ E(X) - +∞ 0 f (x) log 1 -θ + 2 θ e -λ x dx. ( 15
)
Let us now expand the two integrals by using the logarithmic expansion: log(1

-u) = - +∞ m=1 u m m , |u| < 1. Since | γ 1+γ e -λ x (1 -θ + θe -λ x )| < 1, we have +∞ 0 f (x) log 1 - γ 1 + γ e -λ x (1 -θ + θe -λ x ) dx = - +∞ m=1 1 m γ 1 + γ m +∞ 0 f (x) e -λmx (1 -θ + θ e -λ x ) m dx = - +∞ m=1 m =0 1 m m (1 -θ) m- γ 1 + γ m θ +∞ 0 f (x) e -λ( +m)x dx = +∞ m=1 m =0 R m, ,
where

R m, = - 1 m m (1 -θ) m- γ 1 + γ m θ M X [-λ( + m)],
M X (t) denotes the moment generating function defined by [START_REF] Nadarajah | An extension of the exponential distribution[END_REF].

For the second integral in [START_REF] Percontini | The G-negative binomial family: General properties and applications[END_REF], since

|θ(1 -2e -λ x )| < 1, we have +∞ 0 f (x) log 1 -θ + 2θe -λ x dx = - +∞ m=1 θ m m +∞ 0 f (x)(1 -2e -λ x ) m dx = - +∞ m=1 m =0 1 m m θ m (-2) +∞ 0 f (x) e -λ x dx = +∞ m=1 m =0 U m, ,
where

U m, = - 1 m m θ m (-2) M X (-λ ).
Entropy 2. Let us now focus our attention on the Rényi entropy [START_REF] Renyi | On measures of entropy and information[END_REF]:

J R (β) = 1 1-β log +∞ -∞ [f (x)] β dx .
Similar mathematical arguments used for [START_REF] Gupta | Generalized Exponential Distributions[END_REF] give :

1 - γ 1 + γ e -λ x (1 -θ + θe -λ x ) β γ = +∞ i=0 i k=0 C i,k e -λ x (i+k) ,
where

C i,k = β/γ i i k γ θ -γ 1 + γ i θ 1 -θ k .
On the other hand, observing that |θ(1 -2e -λ x )| < 1, similar mathematical arguments used for [START_REF] Gupta | Generalized Exponential Distributions[END_REF] give :

1 -θ + 2 θ e -λ x β = 1 -θ(1 -2 e -λ x ) β = +∞ j=0 j =0 D j, e -λ x ,
where

D j, = β j j (-θ) j (-2) .
Hence [f (x)] β can be expanded as

[f (x)] β = +∞ i=0 i k=0 +∞ j=0 j =0 F i,k,j, e -λ x (i+k+ +β) ,
where

F i,k,j, = λ (1 + γ) (1 + γ) 1 γ +1 -1 β C i,k D j, .
Hence

+∞ -∞ [f (x)] β dx = +∞ 0 +∞ i=0 i k=0 +∞ j=0 j =0 F i,k,j, e -λ x (i+k+ +β) dx = 1 λ +∞ i=0 i k=0 +∞ j=0 j =0 F i,k,j, 1 i + k + + β . Therefore J R (β) = 1 1 -β log +∞ -∞ [f (x)] β dx = 1 1 -β   -log(λ) + log   +∞ i=0 i k=0 +∞ j=0 j =0 F i,k,j, 1 i + k + + β     .
Entropy 3. We now focus our attention on the entropy introduced by [START_REF] Mathai | On generalized distributions and pathways[END_REF]:

J M H (δ) = 1 δ-1 +∞ -∞ [f (x)] 2-δ dx -1 .
Proceeding as for J R (β) with 2 -δ instead of β, we obtain

[f (x)] 2-δ = +∞ i=0 i k=0 +∞ j=0 j =0 G i,k,j, e -λ x (i+k+ +2-δ) ,
where

G i,k,j, = λ (1 + γ) (1 + γ) 1 γ +1 -1 2-δ (2 -δ)/γ i i k γ θ -γ 1 + γ i θ 1 -θ k 2 -δ j j (-θ) j (-2) .
Hence

J M H (δ) = 1 δ -1 +∞ -∞ [f (x)] 2-δ dx -1 = 1 δ -1   1 λ +∞ i=0 i k=0 +∞ j=0 j =0 G i,k,j, 1 i + k + + 2 -δ -1   .
Some numerical values for the three entropies are given in Table 1. It can be observed that these entropies decrease with increasing the parameter values. Moreover, one can see that J M H (δ) has the smallest values comparing with the other entropies considered here.

Conditional moments and mean deviations

Here, we introduce an important lemma which will be used in the next sections.

Lemma 1. Let J r (t) = t 0 x r f (x)dx and Γ(t, ν) = t 0 x ν-1 e -x dx be the lower incomplete Gamma function. Then we have Proof. Using the equation ( 8), we have

J r (t) = 1 λ r+1 +∞ i=0 i k=0 B i,k (1 -θ) Γ{λ(i + k + 1)t, r + 1} (i + k + 1) r+1 + 2 θ Γ{λ(i + k + 2)t, r + 1} (i + k + 2) r+1 . ( 16
)
J r (t) = t 0 x r f (x) dx = +∞ i=0 i k=0 B i,k (1 -θ) t 0
x r e -λx(i+k+1) dx + 2θ t 0

x r e -λx(i+k+2) dx

= 1 λ r+1 +∞ i=0 i k=0 B i,k (1 -θ) Γ{λ(i + k + 1)t, r + 1} (i + k + 1) r+1 + 2 θ Γ{λ(i + k + 2)t, r + 1} (i + k + 2) r+1 .
The r-th conditional moments of the NWTE distribution is given by

E(X r | X > t) = 1 1 -F (t) +∞ t x r f (x) dx = 1 S(t) [E(X r ) -J r (t)] . (17)
It can be expressed using ( 5), ( 13) and Lemma 1. The same remark holds for the r-th reversed moments of the NWTE distribution given by

E(X r | X ≤ t) = 1 F (t) t 0 x r f (x)dx = 1 F (t) J r (t).
The mean deviations of X about the mean µ = E(X) can be expressed as δ = 2µF (µ) -2J 1 (µ) and the mean deviations of X about the median M has the form η = µ -2J 1 (M ).

(Reversed) Residual life functions

Residual lifetime function

The residual life is described by the conditional random variable R (t) = X -t | X > t, t ≥ 0. Using [START_REF] Marshall | A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families[END_REF], the survival function of the residual lifetime R (t) for the NWTE distribution is given by

S R (t) (x) = S(x + t) S(t) = 1 -1 -γ 1+γ e -λ(x+t) 1 -θ + θ e -λ(x+t) 1 γ +1 1 -1 -γ 1+γ e -λt (1 -θ + θ e -λ t ) 1 γ +1 , x ≥ 0.
The associated cdf is given by

F R (t) (x) = 1 -γ 1+γ e -λ(x+t) 1 -θ + θ e -λ(x+t) 1 γ +1 -1 -γ 1+γ e -λt (1 -θ + θe -λt ) 1 γ +1 1 -1 -γ 1+γ e -λt (1 -θ + θe -λt ) 1 γ +1
. The corresponding pdf is given by

f R (t) (x) = λ 1 -γ 1+γ e -λ(x+t) 1 -θ + θ e -λ(x+t) 1 γ e -λ (x+t) 1 -θ + 2 θ e -λ (x+t) 1 -1 -γ 1+γ e -λt (1 -θ + θe -λt ) 1 γ +1
.

The associated hrf is given by

h R (t) (x) = λ 1 -γ 1+γ e -λ(x+t) 1 -θ + θ e -λ(x+t) 1 γ e -λ (x+t) 1 -θ + 2 θ e -λ(x+t) 1 -1 -γ 1+γ e -λ(x+t) 1 -θ + θ e -λ (x+t) 1 γ +1
.

The mean residual life is defined as

K(t) = E(R (t) ) = 1 S(t) +∞ t xf (x) dx -t = 1 S(t) [E(X) -J 1 (t)] -t,
where f (x) is given by (4), S(t) is given by ( 9), E(X) is given by ( 13) and J 1 (t) is given by Lemma 1. Further, the variance residual life is given by

V (t) = V ar(R (t) ) = 2 S(t) +∞ t xS(x)dx -2tK(t) -[K(t)] 2 = 1 S(t) E(X 2 ) -J 2 (t) -t 2 -2tK(t) -[K(t)] 2 ,
where E(X 2 ) is given by( 13) and J 2 (t) is given by Lemma 1. Some numerical values for the mean residual life are displayed in Table 2 for various choices of the parameters γ and θ at the time points t = 1, 3, 5, 7, 10. It can be seen that, the mean residual life increases with increasing the time points t, also decreases with increasing γ and θ.

Reversed residual life function

The reverse residual life is described by the conditional random variable R (t) = t -X | X ≤ t, t ≥ 0. Using (3), the survival function of the reversed residual lifetime R (t) for the NWTE distribution is given by

S R (t) (x) = F (t -x) F (t) = 1 -γ 1+γ e -λ(t-x) 1 -θ + θe -λ(t-x) 1 γ +1 1 -γ 1+γ e -λt (1 -θ + θe -λ t ) 1 γ +1 , 0 ≤ x ≤ t.
The associated cdf is given by

F R (t) (x) = 1 -γ 1+γ e -λt (1 -θ + θe -λt ) 1 γ +1 -1 -γ 1+γ e -λ(t-x) 1 -θ + θ e -λ (t-x) 1 γ +1 1 -γ 1+γ e -λt (1 -θ + θe -λt ) 1 γ +1
. The corresponding pdf is given by

f R (t) (x) = λ 1 -γ 1+γ e -λ(t-x) 1 -θ + θ e -λ (t-x) 1 γ e -λ (t-x) 1 -θ + 2 θ e -λ (t-x) 1 -γ 1+γ e -λ t (1 -θ + θ e -λ t ) 1 γ +1
.

The associated hrf is given by

h R (t) (x) = λ 1 -γ 1+γ e -λ(t-x) 1 -θ + θ e -λ (t-x) 1 γ e -λ(t-x) 1 -θ + 2 θ e -λ (t-x) 1 -γ 1+γ e -λ (t-x) 1 -θ + θ e -λ (t-x) 1 γ +1
.

Moreover, the mean reversed residual life is defined as

L(t) = E(R (t) ) = t - 1 F (t) t 0 xf (x)dx = t - J 1 (t) F (t) ,
where f (x) is given by (4), F (t) is given by (3) and J 1 (t) is given by Lemma 1. Also, the variance reversed residual life is given by

W (t) = V ar(R (t) ) = 2tL(t) -(L(t)) 2 - 2 F (t) t 0 xF (x)dx = 2tL(t) -(L(t)) 2 -t 2 + J 2 (t) F (t) ,
where J 2 (t) is given by Lemma 1.

In Table 3, we give some numerical values for the mean reversed residual life with different choices of the parameters γ and θ at the time points t = 1, 3, 5, 7, 10. From this table, the mean reversed residual life increases with increasing the time points t and with increasing γ and θ.

Estimation

When the parameters λ, γ and θ of the NWTE distribution need to be estimated, several estimation approaches are possible. In this section, we investigate the maximum likelihood estimates (MLEs) of these parameters. Then we propose three goodness-of-fit statistics to compare the densities fitted to any data set.

Maximum likelihood estimation

Let (x 1 , . . . , x n ) be a random samples of size n from the NWTE distribution. Set Θ = (λ, γ, θ) T , then the MLE of Θ can be determined by maximizing the log-likelihood function (Θ) given by

(Θ) = n log(λ) + log(γ + 1) -log [1 + γ] 1 γ +1 -1 + 1 γ n i=1 log 1 + γ -γe -λxi 1 -θ + θe -λxi -λ n i=1 x i + n i=1 log 1 -θ + 2θe -λxi .
Alternatively, by differentiating (Θ), the MLE of Θ can be obtained by solving the nonlinear loglikelihood system equations given by

∂ (Θ) ∂λ = n λ - n i=1 x i - n i=1 2e -λxi θx i 1 -θ + 2e -λxi θ + n i=1 2e -2λxi θ + e -λxi (1 -θ) 1 + γ -e -λxi γ (1 -θ + e -λxi θ) x i = 0, ∂ (Θ) ∂γ = n   1 1 + γ + (1 + γ) 1+ 1 γ {-γ + log(1 + γ)} γ 2 -1 + (1 + γ) 1+ 1 γ   + 1 γ n i=1 1 -e -λxi 1 -θ + e -λxi θ 1 + γ -e -λxi γ (1 -θ + e -λxi θ) - 1 γ 2 n i=1 log 1 + γ -e -λxi γ 1 -θ + e -λxi θ = 0, ∂ (Θ) ∂θ = n i=1 -1 + 2 e -λxi 1 -θ + 2 e -λxi θ - n i=1 e -λxi -1 + e -λxi 1 + γ -e -λxi γ (1 -θ + e -λxi θ) = 0.
By solving the equations above simultaneously, we can obtain the MLE Θ of Θ, with components providing the MLEs λ, γ, θ of λ, γ, θ respectively. Various numerical iterative techniques can be used for estimating these parameters. In this study, we consider the iterative algorithm inherent to the NMaximize command in the symbolic computational package Mathematica.

Under some regularity conditions, the asymptotic normality of the MLEs is guaranteed; the asymptotic distribution of ( Θ -Θ) is N 3 (0 3 , I(Θ) -1 ), where I(Θ) = E(J(Θ)) denotes the expectation of the information matrix: J(Θ) = {J rs (Θ)}, (r, s) ∈ {λ, γ, θ}. Thus confidence intervals or Wald test can be constructed for the parameters.

Goodness-of-Fit Statistics

In order to evaluate the goodness-of-fit of the fitted models, we consider the Anderson-Darling statistics (A * ), the Cramér-von Mises statistics (W * ) and the Kolmogrov-Simnorov statistics (K-S), given by

A * = 2.25 n 2 + 0.75 n + 1 -n - 1 n n i=1 (2i -1) log (z i [1 -z n-i+1 ]) , W * = 0.5 n + 1 n i=1 z i - 2i -1 2n 2 + 1 12n , K-S = Max i n -z i , z i - i -1 n ,
where z i = cdf(y i ) and the y , i s are the ordered observations. The associated P -values are determined. The better distribution in terms of fit is the one having the smallest statistics and largest P -values.

Applications

This section is devoted to the data analyses of two data sets in environmental sciences, namely hydrology, where we compare the fit of our new distributions and some well-known distributions. The best model among them is then selected.

Data fitting

We consider the data sets: "Ground-water data (GWD)" described in Table 1 of Bhaumik and Gibbons [START_REF] Bhaumik | One-Sided Approximate Prediction Intervals for at Least p of m Observations From a Gamma Population at Each of r Locations[END_REF] and "Flood data (FD)" described in (Akinsete et al. [START_REF] Akinsete | The beta-Pareto distribution[END_REF]).The data of GWD represent vinyl chloride concentrations (n = 34) collected from clean upgradient monitoring wells. The data of FD represent flood rates (for the years 1935-1973) (n = 39) for the Floyd River located in James, Iowa, USA. The descriptive statistics of both data sets are summarized in Table 4. From this table, the data are over-dispersed and having skewness and kurtosis. For each data set, the NWTE model is compared with the following distributions.

• The Gamma distribution with pdf given by

f (x) = x k-1 e -λ -1 x λ k Γ(k) , x > 0, k > 0, λ > 0.
• The Marshal-Olkin exponential distribution (MOE) [START_REF] Marshall | A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families[END_REF] with a pdf given by

f (x) = λ β e -λx [1 -(1 -β) e -λ x ] 2 , x > 0 , β > 0 , λ > 0.
• The Nadarajah-Haghighi exponential distribution (NHE) [START_REF] Nadarajah | An extension of the exponential distribution[END_REF] with a pdf given by

f (x) = αλ(1 + λx) α-1 e 1-(1+λx) α , x > 0 , λ > 0, α > 0.
• The exponentiated exponential distribution (EE) [START_REF] Gupta | Generalized Exponential Distributions[END_REF] with a pdf given by f (x) = αλ 1 -e -λx α-1 e -λx , x > 0, λ > 0, α > 0.

• The transmuted Weibull distribution (TW) [START_REF] Aryal | transmuted Weibull distribution: a generalization of the Weibull probability distribution[END_REF] with a pdf given by

f (x) = η σ x σ η-1 e -( x σ ) η 1 -λ + 2 λ e -( x σ ) η , x > 0, η > 0, σ > 0, λ ∈ [-1, 1].
• The transmuted generalized exponential distribution (TGE) [START_REF] Khan | Transmuted generalized exponential distribution: A generalization of the exponential distribution with applications to survival data[END_REF] with a pdf given by

f (x) = α θ e -θ x 1 -e -θ x α-1 1 + λ -2 λ 1 -e -θ x α , x > 0, α > 0, θ > 0, λ ∈ [-1, 1].
• The transmuted linear exponential distribution (TLE) [START_REF] Tiana | Transmuted Linear Exponential Distribution: A New Generalization of the Linear Exponential Distribution[END_REF] with a pdf given by

f (x) = (β + θx) e -(β x+ θ 2 x 2 ) 1 -λ + 2 λ e -(βx+ θ 2 x 2 ) , x > 0, β > 0, θ > 0, λ ∈ [-1, 1].
• The Kappa distribution [START_REF] Mielke | Another Family of Distributions for Describing and Analyzing Precipitation Data[END_REF] with a pdf given by

f (x) = α θ β x β θ-1 α + x β αθ -(α+1) α , x > 0, α > 0, β > 0, θ > 0.
The MLEs with their standard errors are given in Tables 5 and6 for both data sets along with the goodness-of-fit statistics for each distribution. We can see in Tables 5 and6 that the NWTE distribution has the smallest statistics and the largest P -value; it provides the best fit among the considered distributions. This conclusion is confirmed again by Figure 4.

Hydrologic parameters

The nice fit properties of the NWTE distribution motivates the determination of three important hydrologic parameters for the considered data sets: the return level, the conditional mean of the event data and the mean deviation about the return level.

Return level

A return period is an estimate of the likelihood of an event, such as a flood or a river discharge flow to occur. The probability, return period and return level of flood data and ground water contamination data can be estimated using the equation; P (x T ) = 1 -F (x T ), T = 1/P (x T ) and x T = F -1 1 -1 T , respectively, where F -1 (•) is the inverse of the cdf F (x) and P (x T ) called exceedance probability (see, for instance, [START_REF] Ahammed | Variability of annual daily maximum rainfall of Dhaka, Bangladesh[END_REF][START_REF] Nadarajah | Maximum daily rainfall in South Korea[END_REF]). The return level x T under the NWTE distribution is obtained by where x T > 0 and T ≥ 1. Table 7 provides estimates of the return level x T of the ground water contamination data and flood data, respectively, for the return periods T = 2, 5, 10, 20, 50, 100, 200 years based on replacing the parameters λ, γ, θ by their ML estimates in Tables 5 and6. Moreover, the return periods for some largest values of the both data sets are reported in Table 8 and computed using T = 1/P (x T ), where P (x T ) = S(x T ) is the survival function of NWTE distribution given by (1 + γ) (1 + γ)

x T = - 1 λ ln            1 -       1 + θ -(1 + θ) 2 -4 θ γ 1 -1 T [1 + γ] 1 γ +1 -1 + 1 γ γ+1 -1 2 θ                  ,
1 γ +1 -1
, where λ, γ, θ, are the ML estimates corresponding the used data and are given in Tables 5 and6.

Conditional mean of the event data

The conditional mean of the event (GWD or FD) data based on equation ( 17) is defined as

E(X | X > Q) = 1 S(Q) ∞ Q x f (x) dx,
where S(x) is the survival function of NWTE distribution and Q is a value of the event. For example, for the GWD E(X | X > 8.0 m3/s)=10.1378 and FD E(X | X > 71500 mm) = 81788.2.

Mean deviation about the return level

The mean deviation about the return level is the mean of the distances of each value from their return level and it is a measure of the scatter in a population. The mean deviation about return level can be defined as where m(x T ) =

ξ = ∞ 0 |x -x T | f (x)dx = 2x T F (x T ) -x T + µ -2m(x T ),
x T 0 x f (x)dx and f (x) is the density of the NWTE(λ, γ, θ). Table 7 provides mean deviation about the return level, m(x T ), for the return periods T = 2, 5, 10, 20, 50, 100, 200 for the GWD and FD, respectively, noting that we replace the parameters in f (x) (4) by their ML estimates for the corresponding data.

1 Figure 1 :

 11 Figure 1: Plots of the NWTE density.

Figure 2 :

 2 Figure 2: Plots of the NWTE hrf.

Figure 3 :

 3 Figure 3: Plots of the Skewness and Kurtosis of NWTE for λ = 1.

Figure 4 :

 4 Figure 4: Plots of the NWTE density and cdf estimates superimposed on the histogram and empirical cdf for the used data sets.

1 γ

 1 +1 1 -1 -γ 1+γ e -λx 1 -θ + θe -λx 1 γ +1

Table 1 :

 1 Entropy for several arbitrary parameter values with λ = 1

	γ ↓ θ = 0.5	H(f )	J R (0.5)	J M H (0.5)
	0.1	0.94207	1.32902	0.63437
	0.4	0.91244	1.30828	0.61079
	0.8	0.88415	1.28871	0.58774
	1.2	0.86349	1.27456	0.57056
	1.5	0.85122	1.26622	0.56021
	1.8	0.84092	1.25926	0.55144
	2.0	0.83492	1.25522	0.54629
	θ ↓ γ = 0.7
	-0.9	1.40499	1.66797	0.94970
	-0.5	1.33563	1.61265	0.91027
	-0.2	1.24327	1.54839	0.84921
	0.1	1.11731	1.46168	0.75998
	0.4	0.95418	1.34297	0.64014
	0.6	0.82110	1.23549	0.54159
	0.8	0.66351	1.08602	0.42654

Table 2 :

 2 Mean residual life function for arbitrary parameter values with λ = 1

	γ ↓ θ = 0.5	t → 1.0	3.0	5.0	7.0	10
	0.1	0.918095	0.982144	0.997423	0.999648	0.999982
	0.7	0.900668	0.980089	0.997152	0.999612	0.999981
	1.1	0.894316	0.979369	0.997057	0.999599	0.999980
	1.6	0.889012	0.978778	0.996979	0.999588	0.999979
	2.0	0.885996	0.978447	0.996936	0.999582	0.999979
	θ ↓ γ = 0.5	t → 1.0	3.0	5.0	7.0	10
	-0.9	1.219460	1.027861	1.003735	1.000504	1.000025
	-0.5	1.162075	1.020915	1.002810	1.000380	1.000018
	-0.1	1.088250	1.011465	1.001542	1.000208	1.000010
	0.1	1.040826	1.004807	1.000638	1.000086	1.000004
	0.5	0.905030	0.980593	0.997218	0.999621	0.999981
	1.0	0.511566	0.500207	0.500004	0.500000	0.500001

Table 3 :

 3 Mean reversed residual residue life function for arbitrary parameter values with λ = 1

	γ ↓ θ = 0.5	t → 1.0	3.0	5.0	7.0	10
	0.1	0.568855	2.180127	4.080871	6.059352	9.054745
	0.7	0.579174	2.217077	4.126584	6.107094	9.102934
	1.1	0.583962	2.232563	4.145442	6.126726	9.122736
	1.6	0.588572	2.246646	4.162445	6.144397	9.140554
	2.0	0.591508	2.255231	4.172744	6.155087	9.151330
	θ ↓ γ = 0.5	t → 1.0	3.0	5.0	7.0	10
	-0.9	0.415915	1.657881	3.395545	5.331680	8.317248
	-0.5	0.482478	1.812299	3.589379	5.536490	8.524715
	-0.1	0.528625	1.968720	3.792343	5.751628	8.742699
	0.1	0.546765	2.047690	3.897321	5.863170	8.855739
	0.5	0.576254	2.207155	4.114411	6.094402	9.090128
	1.0	0.604164	2.409363	4.399438	6.399130	9.399122

Table 4 :

 4 Descriptive statistics of both data sets

	Data	Mean	Median	SD	Kurtosis	Skewness	M1	M2	Min	Max
	GWD	1.87941	1.15	1.95259	5.00541	1.60369	1.45692	0.8	0.1	8.0
	FD	6771.1	3570	11695.7	25.4436	4.55806	5872.77	2180	318	71500
	SD = Standard Deviation , M1 = Mean deviation about the mean,
	M2 = Mean deviation about the median

Table 5 :

 5 Comparison of fit of the NWTE using different methods of estimation for GWD

				MLE's
	Distributions	Estimates			A *	W *	KS	P -value
	Gamma(k, λ)	1.062685	1.768549		0.320322	0.051617	0.097341	0.904069
		(0.228152)	(0.480351)
	MOE(β, λ)	0.822837	0.481811		0.246723	0.0322959	0.0876375	0.956463
		(0.484526)	(0.173277)
	NHE(α, λ)	0.900308	0.631993		0.24749	0.0332962	0.0838071	0.97073
		(0.344202)	(0.415966)
	EE(α, λ)	1.076412	0.558049		0.325543	0.0528397	0.0977771	0.901191
		(0.247363)	(0.124162)
	TW(η, σ, λ)	1.076390	2.392822	0.418645	0.255586	0.0384952	0.083499	0.971723
		(0.146953)	(0.942422)	(0.606964)
	TGE(α, θ, λ)	1.160258	0.480348	0.395341	0.261748	0.0407411	0.0889903	0.950556
		(0.220835)	(0.217127)	(0.501833)
	TLE(β, θ, λ)	0.404133	0.013541	0.391774	0.248782	0.0336508	0.0825401	0.97467
		(0.267787)	(0.047316)	(0.373851)
	Kappa(α, θ, β)	1.428222	1.236928	1.304859	0.248987	0.037697	0.0871235	0.958587
		(1.0934)	(0.565213)	(0.489635)
	NWTE(λ, γ, θ)	0.465010	9.179478	0.344129	0.234947	0.0320984	0.0793788	0.982912
		(0.204453)	(50.0745)	(0.75164)

Table 6 :

 6 Comparison of fit of the NWTE using different methods of estimation for FD

				MLE's
	Distributions	Estimates			A *	W *	KS	P -value
	Gamma(k, λ)	0.919695	7362.32		1.2662	0.210505	0.147184	0.366821
		(0.182011)	(1906.34)
	MOE(β, λ)	0.293231	0.000069		1.08796	0.148357	0.142366	0.407993
		(0.205071)	(0.000038)
	NHE(α, λ)	0.609712	0.000374		0.900554	0.117556	0.136633	0.460365
		(0.127014)	(0.000163)
	EE(α, λ)	0.968901	0.000144		1.34484	0.23414	0.150306	0.341607
		(0.212611)	(0.000032)
	TW(η, σ, λ)	0.961730	10522.55	0.805980	0.780029	0.102975	0.11102	0.722328
		(0.105471)	(2312.24)	(0.209698)
	TGE(α, θ, λ)	1.081744	0.000103	0.800145	0.795675	0.121992	0.108377	0.749403
		(0.215061)	(0.000032)	(0.21265)
	TLE(β, θ, λ)	0.000095	8.1×10 -12	0.801661	0.782863	0.110629	0.108711	0.74601
		(0.000021)	(1.3×10 -9 )	(0.207988)
	Kappa(α, θ, β)	0.038151	27.732540	1496.464	3.2369	0.658064	0.218663	0.048011
		(0.095815)	(67.9672)	(253.781)
	NWTE(λ, γ, θ)	0.000097	29.109413	0.808284	0.775836	0.111574	0.108037	0.75285
		(0.000022)	(154.575)	(0.199919)

Table 7 :

 7 Return level estimates xT for T and mean deviation about it

		GWD		FD
	T	x T	ξ	x T	ξ
	2	1.24872	1.31433	4143.93	4422.2
	5	2.97016	1.9147	9802.81	6392.9
	10	4.34956	2.89705	14378.1	9653.1
	20	5.7765	4.11844	19306.6	13873.8
	50	7.70555	5.92134	26491.6	20592.4
	100	9.18173	7.35495	32468.3	26397.4
	200	10.665	8.81682	38858.8	32695.9

Table 8 :

 8 Return periods for some largest values of the GWD and FD

	Values of the GWD	Return period	Values of the FD	Return period
	4.0	8.41121	13900	9.32193
	5.1	14.4311	15100	11.1074
	5.3	15.8987	17300	15.1842
	6.8	32.5876	20600	23.7727
	8.0	57.4359	71500	5183.11


Irina Heim
email: i.heim@pgr.reading.ac.uk

Value Co-creation and Local Content Development: Transformation, Digitalization and Innovation in the Oil and Gas Industry

Keywords: Oil and gas industry, Value Co-creation Local content development, Digitalization, Clusters

This study explains how the co-creation of value in networks can lead to technological upgrades in a local industry (local content development) through transformation of business processes, digitalization of drilling rigs and innovation in oil extraction in the oil and gas industry of Kazakhstan. Theoretical perspectives on local content development (LCD) are predominantly informed by economic and political perspectives. The aim of this paper is to develop a strategic perspective on LCD in clusters. This is qualitative research which uses a case study approach.

Introduction

The broad definition of local content policy (LCP) assumes that it is -"an industrial tool that can enable domestic producers to expand their activities, at least partially with domestic inputs, and gain access to international technological and managerial expertise… [in order to] enhance their competitiveness" [START_REF] Hansen | Toward a strategic management perspective on local content in African extractives[END_REF][START_REF] Kalyuzhnova | Local content policies in resource-rich countries[END_REF]. While literature on this topic is currently emerging, this research agenda is predominantly informed by economic and political perspectives, strategic perspectives are virtually absent [START_REF] Hansen | Toward a strategic management perspective on local content in African extractives[END_REF]. This is in itself problematic, because the reason why LCPs may fail is that they are based on an insufficient understanding of stakeholders' strategies and interests. This paper fills this gap in theory by embracing the role of LCPs in technological upgrade. Technological upgrade is defined as increased organizational performance and competitiveness as a result of improved technology capabilities.

2 Collaborative Approach to LCP Shapiro and Rabinowitz [START_REF] Shapiro | Punishment versus cooperation in regulatory enforcement: a case study of OSHA[END_REF] provided an explanation for cooperative approach to regulation in economics, suggesting that collaborative techniques have to be combined with punishment. A collaborative approach to LCPs that defines the expectations of the government, while providing the international oil companies with flexibility to develop its own local content plans and procurement procedures can achieve greater results. The reason why companies collaborate is joint creation of value, or value co-creation. In the management literature, value co-creation is a paradigm that has emerged from the service management field, innovation management studies, and marketing and consumer research [START_REF] Galvagno | Theory of value co-creation: a systematic literature review[END_REF]. It states that organizations interact with each other for the development of new business opportunities. Purposeful interaction creates benefitsdriving dialogue, learning, and resource transfer. Firms act as resource integrators, as specialization forces them to access existing knowledge, skills, competences, people, products, and available investment [START_REF] Gummesson | Marketing as value co-creation through network interaction and resource integration[END_REF]. This interaction which technological platforms often mediate, leads to innovation, participation, and improved services [START_REF] Galvagno | Theory of value co-creation: a systematic literature review[END_REF]. Therefore, ICT technologies has effect on performance of local companies in the oil and gas cluster, a network of interconnected international and local companies, including small and medium enterprises, specialized suppliers, service providers, firms in related industries, associated institutions (universities, standard agencies, and trade associations), government and citizens that co-create value and develop local content through interactions and exchange of resources, technology and management skills with each other. The role of LCP is to support value co-creation in clusters. This can lead to the technological upgrade, i.e. organizational performance and competitiveness based on improved technology capabilities and further local content development.


Yutong Li
email: y.li8@pgr.reading.ac.uk

The Interplay of FDI and R&D: A Study in the Seven Developed Countries

Keywords: FDI, R&D, Developed countries

With the globalization process in the recent two decades, more and more countries are focusing on improving their technology to get more comparative advantages than other countries. The purpose of this paper is trying to explain the causality between foreign direct investment and R&D. We will analysis seven developed countries (Denmark, Finland, Ireland, Japan, the Netherland, the United Kingdom, and the United States) and use pooling data analysis to investigate this causality during the last 35 years.

Introduction

The purpose of this poster paper is to investigate the causality of foreign direct investment and R&D in the developed countries. One of the motivations in the globalisation foreign direct investment, called knowledge, which means the company may separate their headquarters and R&D activities. For example, the headquarters may be in one market, but served the knowledge generated from another market. Therefore, the research question in this paper set as what is the role of research and development played in the foreign direct investment flow? More specifics of this research question could be if one country had high technology, whether that country would like to produce new products in their own country and to attract foreign country investment inward flow; or this country will 'sell' this new technology and conduct investment into other countries?

For example, the United States had the latest technology in the world, which makes it has a more comparative advantage than other countries. Thus, the United States could conduct investment into other countries in some specific area, like agriculture and labour-intensive industry. Therefore, there is a mutual partner relationship between the United States and the other countries. The same situation also could exist in the other nations. We use seven developed countries as an example to estimate if there is a significant relationship between foreign direct investment and local R&D activities.

Pooling Data Analysis

Table 1 indicates a brief regression result in individual countries, which is including three pair relationships between, FDI, local R&D activities, and economic growth. According to the Table, the first causality is about the FDI inward flows and R&D.

In Ireland, FDI inward flows shows a significant negative influence on its R&D. That means in Ireland if they hold a high technology, the most thing they would like do is to conduct foreign investment in the developing country, which is their efficiency seeking motivation, rather than attract foreign investment from another developed country. On the other hand, the R&D in the Netherland displays a positive effect of FDI inward flows, but an antagonistic relationship in the United Kingdom.

For the second causality between FDI outward flows and R&D activities, in Finland, which illustrate a complementary relationship and also they have a positive influence on each other. That means the government prefer to spend more money on their innovation and to improve their unique technology, to have a comparative advantage; rather than conduct foreign investment into other countries. For the last causality, there is only one country, the United States, indicates the economic growth will support the R&D development. However, on the other side of this relationship, local R&D activities show a significant negative influence on GDP growth, which has the same result as the pooling data regression. Again, it has proved that in the short-run, the government could not get benefit from R&D Department.

Conclusion

In this paper, we have analyzed the relationship between foreign direct investment and technology in seven developed countries. In general, according to the pooling data regression (See Table 1), there is a bi-direction relationship between FDI flows and technology. However, the technology has a diminishing changed effect for both FDI inward flows and outward flows, which mean the technology, will not support FDI flows after few years. In another word, because the renewal of technology, the government have to update their technology frequently, to keep the high competition and comparative advantage in the global market.

Table 1 .

 1 Regression result in seven countries

	Country	FDIN to	RD to	FDIO to	RD to	GDP to	RD to
		RD	FDIN	RD	FDIO	RD	GDP
	Denmark	-	+	+ *	-	+	+
	Finland	+	+	+***	+***	-	-
	Ireland	-*	-	-	+	+	+
	Japan	-	-	-	+	-	-
	Netherland	+	+***	+**	-	-	-***
	The United	+	-***	-	-*	-	-***
	Kingdom
	The United	+	-	-**	-***	+*	-***
	States
	* p < 0.1, ** p < 0.05, *** p < 0.01


Emmanuel Bresso

Claire Lacomblez

Anne Pizard

Patrick Rossignol

Faiez Zannad

Malika Smaïl-Tabbone

EB Marie-Dominique Devignes

A data science approach for exploring differential expression profiles of genes in transcriptomic studies-Application to the understanding of ageing in obese and lean rats in the FIGHT-HF project

Keywords: differential expression profiles, heterogeneous graphs, complex network, transcriptomics

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Transcriptomic studies are known to produce huge amounts of information about differentially expressed genes in various situations. The expression levels measured for several thousands of genes in contrasting pairs of situations (different tissues or organs, different physiological state, different age, etc.) allow to calculate foldchange ratios and false discovery rates. Today, data science should be able to derive valuable knowledge units from all these data, by extracting relevant lists of differentially expressed genes and interpreting them. However, in many cases, only a small proportion of all transcriptomic results is finally exploited.

Here we revisit the concept of differential expression profile (DEP) and we propose a combined databasenetwork approach to extract relevant lists of differentially expressed genes based on DEP sharing and to interpret these lists using heterogeneous graph settings before visualisation.

Methods Differential Expression Profile definition.

A transcriptomic study is usually performed on various biological samples derived from tissues or organs under certain conditions. In most cases, expression profiles are defined in a non-supervised manner by clustering the genes on the basis of their expression values across all situations in a study. In other cases, the study involves contrasting situations that lead to the calculation of fold-change (FC) ratio with false discovery rates (FDR) for each gene. Such a setting allows to define differential expression profiles (DEPs).

Let a transcriptomic study be represented by a set of contrasts𝐶 , 𝑖 ∈ {1. . . n}, defined as ordered pairs of situations, and by the FC ratios and FDR values obtained for each gene 𝑔 in each contrast 𝐶 . Four discrete statuses and two modalities have been defined to represent the behaviour of a given gene in a given contrast.

The two modalities correspond to stringent (FC<-2 or FC>2 and FDR <0.05) or not stringent (any FC and FDR <0.05) definitions of differential expression. Under stringent modality, status 𝑆𝑡𝑟1 means differential expression, status 𝑆𝑡𝑟2 and 𝑆𝑡𝑟3 are for up (FC>2) and down (FC<-2) regulation respectively and status 𝑆𝑡𝑟4means no differential expression. Under non stringent modality, similar statuses can be defined and are prefixed by 𝑁𝑆𝑡𝑟.

Formally, a DEP is described by a definition and an extension. The definition of a DEP can be represented as a set of pairs (𝐶 , 𝑠𝑡𝑎𝑡𝑢𝑠) , 𝑖 ∈ {1. . . n} and 𝑠𝑡𝑎𝑡𝑢𝑠 ∈ {𝑆𝑡𝑟1. . . 𝑆𝑡𝑟4, 𝑁𝑆𝑡𝑟1. . . 𝑁𝑆𝑡𝑟4} describing the differential expression status of a gene across a set of contrasts 𝐶 . The extension of a DEP is the set of genes 𝑔 that match the definition. It should be noted that a DEP can also be defined from the differential expression statuses of a given gene across a set of contrasts and then used to select other genes matching the DEP definition.

DEP Relational Data Model.

A relational database named DEPdb was built to store and query transcriptomic results w.r.t. DEPs. The data model covers the contrast definitions between pairs of situations, the transcripts involved and their corresponding genes along with human orthologs when necessary, the FC ratios and associated FDR values, calculated for each transcript and each contrast.

A DEP query interface was implemented on DEPdb to retrieve lists of genes matching a given DEP definition or retrieve a DEP definition for a given gene.

Network-based interpretation of lists of genes.

Various tools exist already to help biologists interpret a list of genes in the light of pathways and interaction networks available in various integrated and curated public sources. Hence, complex graphs are produced showing the multiple interactions existing between the genes of interest and their interactants. Network science should help end-users to interpret such complex networks.

We present here an analysis strategy based on the notion of « heterogeneous graph » in which different types of nodes are interconnected by various types of relationships. Such heterogeneous graph is the basis of the EdgeBox, a « graph knowledge box » constructed by the EdgeLeap company, using the Neo4J graph database system and available public resources [START_REF] Pinet | Integrative network analysis reveals time-dependent molecular events underlying left ventricular remodeling in post-myocardial infarction patients[END_REF]. The EdgeBox currently contains seven types of nodes: protein/gene, disease, pathway, drug, metabolite, gene ontology term and miRNA, and fourteen types of relationships between these nodes (five monopartite, such as protein-protein interactions, and nine bipartite ones). The 2017 version for the FIGHT-HF project concerns human proteins/genes and counts about 211,000 nodes and almost 22 million of relationships.

To interpret a DEP extension, i.e. a list of genes corresponding to a given DEP definition, in the light of the EdgeBox, a first approach consists in querying the EdgeBox for pathways that interconnect at least two genes/proteins of the list. The query result can be completed with the genes/proteins that interact with at least two genes/proteins of the list. When the gene list is short (less than 10 genes), the resulting graph can be interpreted manually and reveals at a glance which pathways or interactions possibly explain why genes share a given DEP definition. This network-query tool has been implemented onto DEPdb connected to the graph knowledge box. It displays the heterogeneous graph using the Cytoscape program for further analyses.

When the gene list grows, the resulting graph becomes extremely complex and impossible to interpret manually.

We propose to filter the nodes for enhancing user interpretation. Enrichment analyses are first carried out on the list of genes to select a small number of nodes corresponding to the top10 of significantly (p<0.001) enriched pathways and biological process GO terms. In parallel, it reveals useful to select subgroups of gene/protein nodes of interest based on their neighbourhood in the EdgeBox. Such selection may be driven by user knowledge. After node reduction, a collection of smaller heterogeneous subgraphs is produced that usually become tractable.

Our combined database-network approach can function iteratively. Starting from a first DEP definition involving a given set of contrasts, the database will return a list of genes that is subsequently displayed as a heterogeneous network. In some cases, the user will select from this network a gene of interest and will query the database to retrieve its DEP across a different set of contrasts. The database will then return on demand all other genes matching this second DEP definition. The new returned extension can then in turn be analysed as a heterogeneous graph, etc.

Results and discussion: Differentially expressed genes upon ageing in heart and kidney of obese and lean SHHF rats Description of the study A transcriptomic study aimed at characterizing simultaneously metabolic syndrome and cardiac, vascular and renal phenotypes in ageing lean and obese SHHF (Spontaneously Hypertensive Heart Failure) rats has been described previously [START_REF] Youcef | Preventive and chronic mineralocorticoid receptor antagonism is highly beneficial in obese SHHF rats[END_REF]. Obesity is induced in SHHF rats by homozygous inactivation of the leptin receptor gene. Rats have been monitored during 11 months (from the age of 1.5 to 12.5 months). In the frame of the FIGHT-HF project, transcriptomic results of this study are newly investigated using our coupled database/network approach.

Definition of two DEPs and extraction of corresponding gene lists

Transcriptomic results have been stored in the DEPdb database as described above. To illustrate our approach we focus on a simple comparison between ageing in obese versus lean rats, in both heart and kidney tissues.

We therefore use four different contrasts from our database: heart samples from « old versus young » lean and obese rats (contrasts 𝐶 and 𝐶 respectively in DEPdb), and kidney samples from « old versus young » lean and obese rats (contrasts 𝐶 and 𝐶 respectively in DEPdb).

We query DEPdb consecutively with two DEP definitions, namely: 𝐷𝐸𝑃 = {(𝐶 , 𝑆𝑡𝑟1), (𝐶 , 𝑆𝑡𝑟4), (𝐶 , 𝑆𝑡𝑟1), (𝐶 , 𝑆𝑡𝑟4)}, and 𝐷𝐸𝑃 = {(𝐶 , 𝑆𝑡𝑟4), (𝐶 , 𝑆𝑡𝑟1), (𝐶 , 𝑆𝑡𝑟4), (𝐶 , 𝑆𝑡𝑟1)} . In other words, 𝐷𝐸𝑃 will retrieve from DEPdb all genes differentially expressed (stringent defintion) upon aging in lean ((𝐶 , 𝑆𝑡𝑟1)and(𝐶 , 𝑆𝑡𝑟1) ) but not in obese ((𝐶 , 𝑆𝑡𝑟4)and (𝐶 , 𝑆𝑡𝑟4)) rats, in heart and kidney samples respectively, and 𝐷𝐸𝑃 will retrieve from DEPdb all genes differentially expressed upon aging in obese ((𝐶 , 𝑆𝑡𝑟1)and (𝐶 , 𝑆𝑡𝑟1)) but not in lean ((𝐶 , 𝑆𝑡𝑟4)and(𝐶 , 𝑆𝑡𝑟4)) rats, in heart and kidney samples respectively. Note that the expression statuses in all other contrasts are not considered here. We retrieved 7 and 55 genes for 𝐷𝐸𝑃 and𝐷𝐸𝑃 definitions, respectively. We also determined with a third appropriate DEP definition that 11 genes are differentially expressed upon aging in both lean and obese rats in heart and kidney samples.

Interpretation with heterogeneous graphs

We subsequently analysed the two extensions of the 𝐷𝐸𝑃 and 𝐷𝐸𝑃 profiles in the light of our network knowledge box. Please note that these gene lists first need to be converted to their human orthologs. The 𝐷𝐸𝑃 extension includes 7 genes and is short enough to be directly analyzed in DEPdb connected with the EdgeBox. We observe that all gene nodes are interconnected together through pathways and GO terms related to cell cycle and cell proliferation, metabolism of proteins and DNA, apoptosis, rhythmic processes, signal transduction and immune system.

Because of its complexity, the𝐷𝐸𝑃 extension of 55 genes underwent the filtering process described above before heterogeneous graph retrieval from the EdgeBox. Interestingly the resulting network clearly reveals two modules of interconnected genes, with 4 genes forming bridges between the two modules. One of these modules involves 22 genes interconnected with the neutrophil degranulation pathway and several GO terms related to inflammation, whereas the other one involves 15 genes interconnected with two pathways : extracellular matrix organization and elastic fiber formation, and with GO terms related to either extracellular matrix organisation or epoxygenase P450 pathway that is involved in anti-inflammatory response.

In summary the comparison of 𝐷𝐸𝑃 and 𝐷𝐸𝑃 extensions reveals that a group of 7 genes involved in particular in cell cycle and cell proliferation is dysregulated in heart and kidney upon ageing in lean rats but not in obese rats. On the contrary, genes that are dysregulated in heart and kidney upon ageing in obese rats but not in lean rats fall into two modules, one related to inflammatory response and the other one related to extracellular matrix organisation. This analysis brings new precise molecular support to the general statement that ageing can have quite different outcomes in lean or obese individuals. The four proteins bridging these two modules (CD44, INTBD2, ANXA2 and CP2E1) can be further investigated to better understand the complementarity between the two groups of biological processes dysregulated upon ageing in the heart and kidney of obese rats. Interestingly, context analysis of the two lists of genes in a more classical framework such as CPDB returned similar results [START_REF] Kamburov | ConsensusPathDB--a database for integrating human functional interaction networks[END_REF].

Conclusion and Perspectives

A formal definition of DEPs has been proposed for transcriptomic studies involving contrasting situations. A coupled database/network approach has been designed and implemented to explore any desired DEP from a study, by extracting the corresponding list of genes and interpreting them using a network knowledge box and heterogeneous graphs. The ageing case study shows that this approach provides useful support for manual interpretation of gene lists and can lead to new hypotheses generation. Obviously, a large number of different DEPs remain to be tested in the same way from the SHHF study.

We are currently searching to design and implement automatic filters to reduce heterogeneous graph complexity in a knowledge-based manner for assisting user interpretation of heterogeneous networks.


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

E Lorin
email: elorin@math.carleton.ca

A simple pseudospectral method for the computation of the time-dependent Dirac equation with Perfectly Matched Layers. Application to quantum relativistic

Keywords: Dirac equation, time-splitting, pseudospectral approximation, Perfectly Matched Layers, high-order accuracy

laser physics.

Introduction

This paper is concerned with the numerical computation of the time-dependent Dirac equation, especially for physical problems involving delocalized wavefunctions. More specifically, we are interested in the numerical solution to the Dirac equation on a truncated domain with absorbing boundary layers, hence limiting the periodic conditions effects. In this goal, we then propose a simple combination of a pseudospectral method with Perfectly Matched Layers (PMLs), allowing to consider delocalized wavefunctions, as for instance observed when quantum relativistic particles are subject to strong fields. Thanks to a relatively new and simple Fourier-based discretization of spatial differential operators, it is possible to impose PML for solving the Dirac equation on a bounded domain. Although, overall, the Fourier-based method still imposes periodic boundary conditions, the outgoing/incoming waves are in fact mainly absorbed. We think that due to the simplicity of the proposed method, most of existing Fourier-based codes could easily be modified to include the proposed methodology. The Dirac equation is a relativistic wave equation which has gain much attention these past 15 years due to the development of 2-d material, such as graphene [START_REF] Fillion-Gourdeau | Numerical computation of dynamical schwinger-like pair production in graphene[END_REF][START_REF] Fillion-Gourdeau | Time-domain quantum interference in graphene[END_REF][START_REF] Katsnelson | Chiral tunnelling and the Klein paradox in graphene[END_REF], intense-laser-molecule interaction [START_REF] Faisal | A four-component Dirac theory of ionization of a Hydrogen molecular ion in a super-intense laser field[END_REF][START_REF] Gonoskov | Probing nonperturbative QED with optimally focused laser pulses[END_REF][START_REF] Mocken | Quantum dynamics of relativistic electrons[END_REF][START_REF] Salamin | Relativistic high-power laser-matter interactions[END_REF], in particular for pair production [START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF], or from heavy ion collisions modeling and simulation for quark-antiquark production [START_REF] Gelis | Quark-antiquark production from classical fields in heavy-ion collisions: 1 + 1 dimensions[END_REF][START_REF] Reinhardt | Theory of positron production in heavy-ion collisions[END_REF][START_REF] Wells | Asymptotic channels and gauge transformations of the time-dependent Dirac equation for extremely relativistic heavy-ion collisions[END_REF]. In the same time, there has been a tremendous progress in the development of efficient and accurate computational methods for the solution to the Dirac equation [START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF][START_REF] Bauke | Accelerating the Fourier split operator method via graphics processing units[END_REF][START_REF] Becker | Solution of the time-dependent Dirac equation by the finite difference method and application for Ca 20+ + U 91+[END_REF][START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF][START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF][START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | Numerical solution of the timeindependent Dirac equation for diatomic molecules: B-splines without spurious states[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF][START_REF] Fillion-Gourdeau | Galerkin method for unsplit 3-d Dirac equation using atomically/kinetically balanced B-spline basis[END_REF][START_REF] Grant | B-spline methods for radial Dirac equations[END_REF][START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF][START_REF] Mocken | Quantum dynamics of relativistic electrons[END_REF][START_REF] Selstø | Solution of the Dirac equation for hydrogenlike systems exposed to intense electromagnetic pulses[END_REF], modeling in particular molecules subject to intense external laser fields. In this spirit, we have developed [START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF], a Schwinger-like pair production procedure from the interaction of intense laser fields with heavy molecules occuring at specific resonances [START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF][START_REF] Fillion-Gourdeau | Enhanced schwinger pair production in many-centre systems[END_REF].

In this paper, we are more specificially interested in the interaction of relativistic interaction of atoms, molecules or wavepackets with intense and short laser pulses. The key point is that the laser field actually delocalizes the wavefunction, and the latter can actually interact with the domain boundary. In order to avoid artificial reflection it is necessary to impose absorbing boundary conditions [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF], absorbing complex potentials or perfectly matched layers [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF]. Theoretically, this approach allows to benefit from the spectral convergence and simplicity of Fourier-based methods, more specifically pseudospectral, on bounded domains reducing the periodic boundary condition effect thanks to artificial wave absorption at the domain boundary. Perfectly Matched Layers are now widely used in many engineering and physics simulations codes [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schroedinger 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] 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] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF][START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[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] Nissenand | An optimized perfectly matched layer for the Schroedinger equations[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][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] to model exterior domains and to avoid any unphysical reflection at the boundary. PML for the Dirac equation were developed in [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF] and approximated with a real space method. The derivation of high-order ABCs for the Dirac equation were proposed in [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF]. We also refer to [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] for an overview of PMLs and ABCs for quantum wave equations including the Dirac equation and to [START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF][START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | An efficient and stable numerical method for the Maxwell-Dirac system[END_REF][START_REF] Bauke | Accelerating the Fourier split operator method via graphics processing units[END_REF][START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF][START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF][START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF][START_REF] Fischer | A B-spline Galerkin method for the Dirac equation[END_REF][START_REF] Guo | Spectral and pseudospectral approximations using Hermite functions: Application to the Dirac equation[END_REF][START_REF] Hammer | A dispersion and norm preserving finite difference scheme with transparent boundary conditions for the Dirac equation in (1+1)d[END_REF][START_REF] Kullie | Dirac-Fock finite element method (FEM) calculations for some diatomic molecules[END_REF][START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF][START_REF] Mocken | FFT-split-operator code for solving the Dirac equation in 2+1 dimensions[END_REF] for different approaches for solving the Dirac equation in real or Fourier space.

The combination of pseudospectral methods and PMLs is possible thanks to the following simple and original idea. For example in the x-direction, let us denote by ξ x the dual Fourier variable and by F x (resp. F -1

x ) the Fourier (resp. inverse Fourier) transform in x. We consider a as a given x-dependent function. Then, for any function f , it is possible to formally rewrite a(x)∂ x f (x) as F -1

x a(x)iξ x F(f )(ξ x ) (x) on an unbounded domain (pseudodifferential operator representation [START_REF] Taylor | Pseudodifferential Operators[END_REF]). Let us remark that the latter is a real space function, although the derivative is approximated using the Fourier transform, and the function (x, ξ x ) → ia(x)ξ x is nothing but the symbol of a(x)∂ x . In practice the function a involves the stretching coordinates function modeling the PML on a bounded domain allowing real space non-reflecting conditions at the domain boundary.

The time-dependent Dirac equation under consideration reads [START_REF] Itzykson | Quantum Field Theory[END_REF]

i∂ t ψ(t, x) = Hψ(t, x), (1)
where ψ(t, x) is the time and coordinate dependent four-spinor, and H is the Hamiltonian operator. The latter is given by

H = α • [cp -eA(t, x)] + βmc 2 + I 4 V (t, x), (2)
where the momentum operator is p = -i∇. More specifically, the Dirac equation that we consider reads [START_REF] Itzykson | Quantum Field Theory[END_REF] i∂ t ψ(t, x) = α x -ic∂ x -eA x (t, x) + α y -ic∂ y -eA y (t, x)

+α z -ic∂ z -eA z (t, x) + βmc 2 + I 4 V (t, x) ψ(t, x), (3)
where ψ(t, x) ∈ L 2 (R 3 ) ⊗ C 4 is the time and coordinate (x = (x, y, z)) dependent fourspinor. In (3), A(t, x) represents the three space components of the electromagnetic vector potential, V (t, x) = eA 0 (t, x) + V nuc. (x) is the sum of the scalar and interaction potentials, e is the electric charge (with e = -|e| for an electron), I 4 is the 4 × 4 unit matrix and α = (α ν ) ν=1,••• ,4 , β are the Dirac matrices. In this work, the Dirac representation is used, where

α ν = 0 σ ν σ ν 0 , β = I 2 0 0 -I 2 . ( 4
)
The σ ν are the usual 2 × 2 Pauli matrices defined as

σ x = 0 1 1 0 , σ y = 0 -i i 0 and σ z = 1 0 0 -1 , (5)
while I 2 is the 2 × 2 unit matrix. Note that the light velocity c and fermion mass m are kept explicit in Eq. ( 2), allowing to adapt the method easily to natural or atomic units (a.u.). The paper is organized as follows. We recall the basics of PMLs for the Dirac equation in Section 2. Section 3 is dedicated to the derivation of the pseudospectral approximation for the Initial Boundary Value Problem (IBVP) time-dependent Dirac equation on a bounded domain with PML. Numerical experiments are presented in Section 4. We finally conclude in Section 5.

PML for the Dirac equation

The time-dependent Dirac equation is then considered on a bounded (truncated) physical domain denoted by D Phy . We add a layer, that is called D PML , surrounding D Phy , and stretch the coordinates in all the directions. The overall computational domain is next defined by: D = D Phy ∪ D PML . We refer to [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] for the construction of PMLs for quantum wave equations and more specifically to [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF] for the derivation and analysis of PMLs for the Dirac equation.

The starting point is the stretching of the real coordinates in the complex plane. In [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF], the author uses the following variables

ν = ν + i ω ν L * ν σ(s)ds,
where ω is the dual Fourier variable to t and the so-called absorbing function σ is such that σ(ν) = 0 if s < L * ν . Then, the partial derivative ∂ ν , with ν = x, y, z, is shown to be formally transformed into

∂ ν ψ(t, •) → ∂ t ∂ t + σ ν (ν) ∂ ν ψ(t, •) = 1 2π R ωe iωt ω -i σ ν (ν) ∂ ν F t ψ(ω, •)dω, (6)
where F t (resp. F ν ) denotes the Fourier transform with respect to t (resp. ν) and

σ ν (ν) = σ(|ν| -L ν ), L * ν |ν| < L ν , 0, |ν| < L * ν .
Let us remark that the transformation (6) can also be formally rewritten as

∂ ν ψ(t, •) → 1 1 + σ ν (ν)∂ -1 t ∂ ν ψ(t, •).
In the following, we will rather consider the following change of variables [START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] involving only the space variable

ν = ν + e iθ ν L * ν σ(s)ds,
where θ ∈ (0, π/2). We then define S ν (ν) := 1 + e iθν σ(ν), with ν = x, y, z. In the following, the partial derivatives are transformed into

∂ ν → 1 S ν (ν) ∂ ν = 1 1 + e iθν σ ν (ν) ∂ ν , (7)
where σ vanishes in D Phy and S ν is equal to 1. From now on, let us consider the transformation [START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF], and the associated new Hamiltonian

H PML = α • [cT • p -eA(t, x)] + βmc 2 + I 4 V (t, x), (8)
where T = (S -1 x (x), S -1 y (y), S -1 z (z)) T . Several types of functions can be selected. An exhaustive study of the absorbing functions σ is proposed in [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF] for Schrödinger equations Type I:

σ 0 (ν + δ ν ) 2 , Type II: σ 0 (ν + δ ν ) 3 , Type III: -σ 0 /ν, Type IV: σ 0 /ν 2 , Type V: -σ 0 /ν -σ 0 /δ ν , Type VI: σ 0 /ν 2 -σ 0 /δ 2 ν ,
for ν = x, y, z. The main difficulty from the pseudospectral point of view is the spacedependence of the coefficients S -1 ν which prevents the direct application of the Fourier transform on the equation. A simple trick will however allow for combining the efficiency of the pseudospectral method and the computation of the non-constant coefficient Dirac equation.

PSeudospectral-PML (PS-PML) method

A splitting of the Dirac Hamiltonian is used below. Although operator splitting is not fundamentally required in the proposed methodology, it allows to simplify the implementation of the method while keeping a good accuracy. Based on (2), we define the operators

A = -icS -1 x (x)α x ∂ x , (9)
B = -icS -1 y (y)α y ∂ y , (10)
C = -icS -1 z (z)α z ∂ z , (11)
D = βmc 2 + I 4 V (t, x) -eα • A(t, x). (12)
From time t n to t n+1 , we then successively solve [START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF] (the x-dependence in the wavefunction argument for notational convenience)

i∂ t ψ (1) (t) = Aψ (1) (t), ψ (1) (t n ) = ψ n , t ∈ [t n , t n 1 ) ( 13
)
i∂ t ψ (2) (t) = Bψ (2) (t), ψ (2) (t n ) = ψ (1) (t n 1 ), t ∈ [t n , t n 2 ) ( 14
)
i∂ t ψ (3) (t) = Cψ (3) (t), ψ (3) (t n ) = ψ (2) (t n 2 ), t ∈ [t n , t n 3 ) ( 15
)
i∂ t ψ (4) (t) = Dψ (4) (t), ψ (4) (t n ) = ψ (3) (t n 3 ), t ∈ [t n , t n+1 ) (16) and ψ n+1 = ψ (4) (t n+1 ), (17)
where

t n i -t n = ∆t for i ∈ {1, • • • , 3}.
In the next subsections, we detail the space-time approximation of equations ( 13)- [START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF]. Cylindrical coordinates could be used as well (see e.g. [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF]).

Pseudospectral approximation in space

We consider for convenience the 3-dimensional system in cartesian coordinates, for ν = x, y, z in the domain

D = [-a x , a x , ] × [-a y , a y ] × [-a z , a z ], i∂ t ψ(t) = -icα ν S -1 ν (ν)∂ ν ψ(t), ψ(t n ) = ψ n , t ∈ [t n , t n+1 )
and we use the same notation as in [START_REF] Antoine | High-order IMEX-spectral schemes for computing the dynamics of systems of nonlinear Schrödinger/GrossPitaevskii equations[END_REF]. We first diagonalize α ν = Π ν D ν Π T ν , where

D ν =     1 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 -1     .
The matrices Π ν are defined as follows

Π x = 1 √ 2     0 1 1 0 1 0 0 -1 1 0 0 1 0 1 -1 0     , Π y = 1 √ 2     0 -i -i 0 1 0 0 1 -i 0 0 i 0 1 -1 0     , Π z = 1 √ 2     1 0 0 -1 0 -1 -1 0 1 0 0 1 0 1 -1 0     .
We set φ := Π T ν ψ, which then satisfies

i∂ t φ(t) = -icD ν S -1 ν (ν)∂ ν φ(t), φ(t n ) = Π T ν ψ n , t ∈ [t n , t n+1 ).
We denote the set of grid-points by

D Nx,Ny,Nz = x k 1 ,k 2 ,k 3 = (x k 1 , y k 2 , z k 3 ) (k 1 ,k 2 ,k 3 )∈O NxNy Nz with O NxNyNz = (k 1 , k 2 , k 3 ) ∈ N 3 / k 1 = 0, • • • , N x -1; k 2 = 0, • • • , N y -1; k 3 = 0, • • • , N z -1 .
Then, let us introduce the following mesh sizes

x k 1 +1 -x k 1 = h x = 2a x /N x , y k 2 +1 -y k 2 = h y = 2a y /N y , z k 3 +1 -z k 3 = h z = 2a z /N z . (18)
The corresponding discrete wavenumbers are defined by ξ := (ξ p , ξ q , ξ r ), where

ξ p = pπ/a x with p ∈ {-N x /2, • • • , N x /2 -1}, ξ q = qπ/a y with q ∈ {-N y /2, • • • , N y /2 -1} and ξ r = rπ/a z with r ∈ {-N z /2, • • • , N z /2 -1}.
In the sequel of the paper, we denote by φ ( ) , with ∈ {1, 2, 3, 4}, the th component of the spectral approximation of φ = Π T ν ψ, with ν = x, y, z. We also use the notation φ ( )

k 1 (t, y, z) = φ ( ) (t, x k 1 , y, z), φ ( ) k 2 (t, x, z) = φ ( ) (t, x, y k 2 , z) and φ ( ) k 3 (t, x, y) = φ ( ) (t, x, y, z k 3 ). The partial Fourier coefficients are such that                        φ ( ) p (t, y, z) = Nx-1 k 1 =0 φ ( ) k 1 (t, y, z)e -iξp(x k 1 +ax) , φ ( ) q (t, x, z) = Ny-1 k 2 =0 φ ( ) k 2 (t, x, z)e -iξq(y k 2 +ay) , φ ( ) r (t, x, y) = Nz-1 k 3 =0 φ ( ) k 3 (t, x, y)e -iξr(z k 3 +az) .
We can then introduce the inverse partial Fourier pseudospectral approximations, in the x-, y-and z-directions, respectively,

                         φ ( ) k 1 (t, y, z) = 1 N x Nx/2-1 p=-Nx/2 φ ( ) p (t, y, z)e iξp(x k 1 +ax) , φ ( ) k 2 (t, x, z) = 1 N y Ny/2-1 q=-Ny/2
φ ( ) q (t, x, z)e iξq(y k 2 +ay) , φ ( )

k 3 (t, x, y) = 1 N z Nz/2-1 r=-Nz/2
φ ( ) r (t, x, y)e iξr(z k 3 +az) .

We finally define the approximate first-order partial derivatives by

                         ∂ x φ ( ) (t n , x k 1 ,k 2 ,k 3 ) ≈ [[∂ x ]]φ ( ) k 1 ,k 2 ,k 3 := 1 N x Nx/2-1 p=-Nx/2 iξ p φ ( ) k 2 ,k 3 p e iξp(x k 1 +ax) , ∂ y φ ( ) (t n , x k 1 ,k 2 ,k 3 ) ≈ [[∂ y ]]φ ( ) k 1 ,k 2 ,k 3 := 1 N y Ny/2-1 q=-Ny/2 iξ q φ ( ) k 1 ,k 3 q e iξq(x k 2 +ay) , ∂ z φ ( ) (t n , x k 1 ,k 2 ,k 3 ) ≈ [[∂ z ]]φ ( ) k 1 ,k 2 ,k 3 := 1 N z Nz/2-1 r=-Nz/2 iξ r φ ( ) k 1 ,k 2 r e iξr(x k 3 +az) . (19)
In the following, the index h will be used (e.g. in

φ n h = {φ n k 1 ,k 2 ,k 3 } k 1 ,k 2 ,k 3 )
to denote a spectral approximation to a given wavefunction (e.g. φ n ). This discretization not only allows to select the spatial step as large as wanted, but it also preserves the very high spatial accuracy, the parallel computing structure and the scalability of the split method developed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF].

Time-Splitting PSeudospectral method with PML (TSSP-PML)

As proposed before, we use a splitting of the TDDE into four time-dependent systems (9), [START_REF] Becker | Solution of the time-dependent Dirac equation by the finite difference method and application for Ca 20+ + U 91+[END_REF], [START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF] and [START_REF] Bermúdez | An exact bounded perfectly matched layer for time-harmonic scattering problems[END_REF]. We denote by x h = (x h , y h , z h ) the nodes of a real-space grid. At iteration n, the approximate wave function is denoted by ψ n h . The TSSP (Time-Splitting pseudoSPectral) algorithm then reads 1. First step: integration of the generalized transport equation in Fourier space in the x-direction. One sets φ n h := Π T x ψ n h , and the system

∂ t φ + cS -1 x (x)Λ x ∂ x φ = 0, φ(t n , •) = Π T x ψ(t n , •),
is approximately solved by

φ n 1 h = φ n h -c∆tS -1 x (x h )Λ x [[∂ x ]]φ n h ,
where [[•]] is defined in [START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[END_REF]. We then deduce ψ n 1 h = Π x φ n 1 h . This gives an approximation of Eq. ( 13), where the corresponding operator is denoted by P (x) h (∆t).

2. Second step: integration of the generalized transport equation in Fourier space in the y-direction

∂ t φ + cS -1 y (y)Λ y ∂ y φ = 0.
We set χ n 1 h := Π T y ψ n 1 h and we approximately solve ( 14) by using the following scheme

φ n 2 h = χ n 1 h -c∆tS -1 y (y h )Λ y [[∂ y ]]χ n 1 h .
We then obtain ψ n 2 h = Π y φ n 2 h . 3. Third step: setting χ n 2 h := Π T z ψ n 2 h , we apply the same procedure in the z-direction, and compute ψ n 3 h . The corresponding operator is denoted by P (z) h (∆t). 4. Fourth step: integration of the source from t n to t n+1

ψ n+1 h = T exp -i t n+1 tn dτ βmc 2 -eα • A h (τ ) × exp -ie t n+1 tn dτ V h (τ ) ψ n 3 h , (20)
where A h is the approximate electric potential. The corresponding operator is denoted by Q n h (∆t). This splitting scheme is fully explicit, and can be compactly expressed as a composition of the four operators

ψ n+1 h = P (x) h (∆t)P (y) h (∆t)P (z) h (∆t)Q n h (∆t)ψ n h .
In practice, a second-order splitting scheme should be implemented, i.e.

ψ n+1 h = P (x) h ∆t 4 P (y) h ∆t 2 P (x) h ∆t 4 P (z) h ∆t 2 Q n h (∆t) P (z) h ∆t 2 P (x) h ∆t 4 P (y) h ∆t 2 P (x) h ∆t 4 ψ n h . (21)
Thanks to the operator [[∂ ν ]], it is possible to implement the real-space approximation of the overall wave function on a bounded domain imposing in particular PMLs, resulting in the TSSP-PML scheme.

Numerical simulation

We propose some numerical experiments illustrating the efficiency of the combined method PML-TSSP method with a second-order operator splitting. We are mainly interested in the evolution of wavepacket either subject to a static potential or of an external field. In our simulations, we have used several libraries which are listed below.

Technical details about the code and parallel computing aspects

The Fourier transforms are performed with the sequential and parallel version 3.3.4 of fftw. Standard linear algebra libraries gsl (version 1.9) and blas were also used. Finally openMPI 1.6.3 was used for message passing with non-blocking communication. The code was implemented in C++ and the compiler is the version 4.7.0 of gcc. In the following, the first three one-dimensional tests are realized by using matlab. The tests are performed by using a C++-code with MPI-library. The method is implemented on the cluser mammouth-parallel II from the RQCHP. The total processing power of this machine is 333 400 GFlops, and possesses 39648 cores: 3216 processors AMD Opteron 12 cores at 2.1 GHz, and 88 processors AMD Opteron 12 cores 2.2 GHz. The total memory is 57.6 TB and the computer-networking communications is Infiniband QDR (4 GB/sec). Regarding the communication between processes, by default mammouth II will select the closest processors in the same node, then the processors from the nodes in the same topological ring, and finally through different switch-levels. Notice that for a low number of processors, the communication between the processes within the same node is done thanks to a shared memory. The detailed of the Dirac equation solver without PML is presented in detailed in [START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF], Section 4. Let us recall the main features of the code in 2-d. In the (x, z) coordinates, we proceed by alternating the directions. We first decompose the domain by layers in the z-direction, then

• We successively perform the evolution (FFT) sequentially in the x-direction, and by layer in the z-direction. Each processor manages one layer in z. A perfect scaling for this step is expected as it does not require any transmission from one node to another.

• For all x, we perform the evolution in the z-direction using the parallel FFT (fftw). The performance of this step is then fully dependent on the parallelization of the one-dimensional FFT.

Notice that the presence of the PML does not deteriorate the efficiency of the overall method compared to usual FFT-methods. In particular, the computational complexity is simple to established, and is given by C after N T time iterations

C = O N T f N x N z log(N x N z ) .
The scalability and performance of the TSSP method without PML is fully studied in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], where in particular it is shown a good speed-up of this method.

Numerical experiments

The second-order TSSP-PML method ( 21) is implemented, where the discrete operators [[•]] given by [START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[END_REF] are used in the directions ν = x, z. We then choose S ν as follows, for ν = x, z,

S ν (ν) = 1, |ν| < L * ν , 1 + e iθ σ(|ν| -L v ), L * ν |ν| < L ν . ( 22
)
Test 1 : wavepacket subject to an external laser field. In this first test, we consider a wavepacket

ψ(x, y, z, 0) = N [1, 0, 0, 0] T × exp -(x 2 + y 2 + z 2 )/∆ 2 exp(ik x x), (23)
subject to an external potential and propagating with a fixed wavenumber k = (5, 0), and ∆ = 128 in [START_REF] Fillion-Gourdeau | Numerical solution of the timeindependent Dirac equation for diatomic molecules: B-splines without spurious states[END_REF]. The time step is given by ∆t = 4.56 × 10 -4 . The physical domain is D = [-8, 8] 2 and is discretized with N x × N y = 256 2 grid points. We impose a linearly polarized electric field such that A x is identically null, and A z (t, x, z) is defined as follows

A z (t, x, z) = A 0 cos(ωx) exp -(T f -t) 2 /2 ,
with A 0 = 100/137, ω = 100, a total of 20 cycles and T f = 1.824. As an example, we report in Fig. 1 the electric potential A z at time T = 0.456 and 1.824 and the initial data. In Fig. 2, we compare the evolution of the density d(T, x, z) = 4 i=1 |ψ i (T, x, z)| 2 as a function of time, with i) periodic boundary conditions, with PML of ii) type III and iii) type IV, with σ 0 = 10 -2 , θ = π/4 as well as iv) type V and VI with δ ν = L ν -L * ν . The snapshots correspond to times T = 0.456, 1.14, 1.368 and 1.824, for respectively the time iterations 1000, 2500, 3000 and 4000. The best absorption is obtained with the Type-VI absorbing function, although the results are relatively close with all the absorbing functions. We report in Fig. 3, the maximum of the density (in logscale) as a function of the time iterations for 6 different absorbing functions, showing that they all provide relatively good properties for truncating the computational domain.

Test 2 : 3-nucleus system. In this test, a Gaussian wavepacket centered at (0, 0) Fig. 4 (Left), is injected in a 3-nucleus interaction potential Fig 4

(Middle) V (x, z) = - Z A (x -x A ) 2 + (z -z A ) 2 + 2 - Z B (x -x B ) 2 + (z -z B ) 2 + 2 - Z C (x -x C ) 2 + (z -z C ) 2 + 2 (24)
with

Z A = Z B = Z C = 10a.u., (x A , z A ) = (-1, 0), (x B , z B ) = (1, 0) and (x C , z C ) = (0, -1).
The overall computational domain is D = [-6, 6] 2 , with N x × N y = 128 2 grid points and ∆t = 3 × 10 -3 . We report on Fig. 4 (Right) the total density at time T = 3 × 10 -1 a.u., i.e. d(T, x, z), for (x, z) ∈ D, where d(T, x, z) = 4 i=1 |ψ i (T, x, z)| 2 , setting ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) as the Dirac wave function. Without external excitation the wavepacket is then trapped by the potential. We then plug an external laser field which drives the wavepacket: a linearly polarized electric field is imposed such that A x is identically null, and A z (t, x, z) is defined as follows

A z (t, x, z) = A 0 cos(ωx) exp -(T -t) 2 /2 ,
with A 0 = 1000/137, ω = 1, and a total of 2 cycles. We report the z-component of the external laser potential in Fig. 5 at final time. We also include a non-null wavenumber in the initial wavepacket, k = (10, 0), out of the ions influence. The wavepacket will eventually be absorbed by the PML unlike the solution without PML.

At times t = 3 × 10 -2 , t = 6 × 10 -2 a.u., t = 3 × 10 -1 and t = 3.6 × 10 -1 , we compare 3 solutions: 1) the solution of reference computed on [-18, 18] × [-18, 18], 2) the solution without PML, and 3) the solution with Type V-PML and σ 0 = 10 -2 , θ = π/4, and δ = 1.2. At time t = 6 × 10 -2 , the 3 solutions look identical, as the wavefunction has not reached the boundary of the computational domain [-6, 6] × [-6, 6]. At t = T = 3 × 10 -1 , we however clearly see that the wavefunction of reference has mainly leaved the region [-6, 6] × [-6, 6], while the PML solution was almost totally absorbed and the solution without PML, was maintained in the computational domain due to non-absobing boundary layers. This is also illustrated on Fig. 7, where the 2 -norm as function of time in the D is represented for the 3 solutions. We see that, unlike the solution without PML, the PML-solution makes decrease the 2 -norm, thanks to the absorbing layers. We also report in Fig. 7 (Left), we represent the solution of reference at time T = 0.3a.u. in a bigger domain [-18, 18] × [-18 , 18]. By comparison, we see that on this example, the solution of reference has also a decreasing This test illustrates that artificial wave-reflecting occurs with real space methods (here a quantum lattice Boltzmann method [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF]) are implemented without PML or ABC. Moreover in Fig. 14 (Right), we observe an exponential convergence up to a certain precision.

Notice that what is reported is not the comparison with a solution of reference, but the maximum of the solution over the physical domain, including the PML. In other words the smaller the maximum, the better the absorption. The saturation error comes from the limit of the PML accuracy and error in time (second-order splitting). This again shows the high accuracy of the TSSP with PML.

Conclusion

In this paper, we have proposed a simple time-splitting pseudospectral method which allows for the numerical computation of IBVP for the Dirac equation with non-reflecting PML. Some numerical experiments have shown that the method preserves the accuracy of the PML with static or time-dependent external potentials. Different absorbing functions were considered and tested. The implementation strategy is developed for the TSSP approximation, but can be directly extended to other numerical schemes including implicit. We think that Fourier-based codes solving the time-dependent Dirac equation can easily be adapted to include the absorbing layers, hence drastically reducing the negative effect of periodic conditions and simultaneously avoiding artificial wave reflections. We next plan to analyze mathematically the stability and convergence of the method as well as the PML accuracy for different types of absorbing functions.

Figure 1 :

 1 Figure 1: Test 1 : initial density (top), A z (T, x, y) at times T = 0.456 (Bottom left) and 1.824 (Bottom right).

2 -

 2 norm in the zone D ( t, d Ref |D (t, •) 2 ), see Fig.7where the 2 -norm is represented in logscale as a function of time (Right). We also report the solution of reference in the domain[-18.18] × [-18, 18] at time T = 0.36a.u.Test 3 : evolution of wavepacket subject to a repulsive weakly nonlinear poten-

Figure 2 :Figure 3 :

 23 Figure 2: Test 1 : density at times T = 0.456, 1.14, 1.368 and 1.824 (first line) with periodic boundary conditions (without PML). The second (respectively third) line reports the results for the density with TSSP-PML, with a PML of Type-III (respectively Type-IV), setting σ 0 = 10 -2 , θ = π/4. The fourth (respectively fifth) line gives the same kind of results but for the Type-V (respectively Type-VI) PML, with σ 0 = 10 -2 , θ = π/4 and δ ν = L ν -L * ν .

Figure 4 :

 4 Figure 4: Test 2 : (Left) Initial density. (Middle) 3-nucleus potential. (Right) Density at time T = 3 × 10 -1 a.u. without external laser field.

Figure 5 :

 5 Figure 5: Test 2 : z-component of the external laser potential: A z (T, x, z) for T = 3 × 10 -1 a.u.

Figure 6 :

 6 Figure 6: Test 2 : From top to bottom: No PML, reference and Type V-PML solutions at time t = 3 × 10 -2 a.u. t = 6 × 10 -2 a.u., t = 3 × 10 -1 a.u. and t = 3.6 × 10 -1 a.u.,

Figure 7 :

 7 Figure 7: Test 2 : (Left) Solution of reference at time T = 6 × 10 -2 a.u. in the domain [-18, 18] × [-18, 18]. (Right) 2 -norm in logscale of the 3 solutions (No PML, Reference, Type V-PML) as a function of time in D.

Figure 8 :

 8 Figure 8: Test 3 : initial density in logscale (Left) and static potential (Right).

Figure 9 :

 9 Figure 9: Test 3 : logarithm of the density at time T = 0.912 (from top-left to bottom-right) for the PML with parameters σ 0 = 1, σ 0 = 10 -1 , σ 0 = 10 -2 , σ 0 = 10 -3 , σ 0 = 5, and finally for a periodic boundary condition (S = 1).

Figure 10 :

 10 Figure 10: Test 3 : logarithm of the density at time T = 0.912 with null Dirichlet boundary condition.This test illustrates that artificial wave-reflecting occurs with real space methods (here a quantum lattice Boltzmann method[START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF]) are implemented without PML or ABC.

Figure 11 :

 11 Figure 11: Test 3 : PML of type IV (Left) maximum of the density as a function of the time iteration with 5 different values of σ and for θ = π/4; (Right) maximum of the density as a function of the time iteration for θ = π/16, π/8, π/4, π/2 and σ 0 = 10 -2 .

Figure 12 :

 12 Figure 12: Test 3 : logarithm of the density at time T = 0.912 for σ 0 = 10 -2 and θ = π/16, θ = π/8, θ = π/4 and θ = π/2.

Figure 13 :

 13 Figure 13: Test 3 : logarithm of the density at time T = 0.912 for σ 0 = 10 -2 , θ = π/4, δ ν = L ν -L * ν , and absorption function σ of types I to IV.

Figure 14 :

 14 Figure 14: Test 3 : initial data (Left), periodic boundary conditions (Middle), convergence with the Type-III PML (Right), for σ 0 = 10 -2 , θ = π/4 and N x = N y = 32, 256, 320, 512, 640, 768.

Acknowledgments. X. ANTOINE was partially supported by the French National Research Agency project NABUCO, grant ANR-17-CE40-0025. E. LORIN received support from NSERC through the Discovery Grant program.

tial. We consider the two-dimensional Dirac equation with a nonlinear potential

with Z A = Z B = 2, (x A , z A ) = (-1, 4) and (x B , z B ) = [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schroedinger equation[END_REF][START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF]. Notice that the sign of the potential makes it repulsive. The initial density with wave vector k = (2, 10) is represented in Fig. 8. More specifically, the initial data is given by

where ∆ = 128 in [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF] and

The overall computational domain is D = [-8, 8] 2 , with N x × N y = 256 2 grid points and ∆t = 4.56 × 10 -4 . We report on Fig. 9 the total density at time T = 0.912 a.u. in logscale, i.e. log(d(T, x, z), for (x, z) ∈ D, where d(T, x, z) =

) as the Dirac wave function. The nonlinearity is numerically treated explicitly. The use of the logscale allows to fairly report the accuracy of the PML. In [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], we take L ν = 8, L * ν = 0.8L ν , θ = π/4 (similarly to the Schrödinger equation [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF]), and we pick the Type-IV PML, i.e. σ : ν → σ 0 /ν 2 (singular profile [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF]). We represent the solution with different values of σ 0 . These results show that the PML properly absorbs except for σ 0 = 5, with reflection magnitude as weak as 10 -6 for σ 0 = 10 -2 . Let us remark that other PMLs could be considered as well [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF], but we intend here to prove the feasibility of the method for various PMLs. For completeness and by comparisons with results obtained in Fig. 9, we also report in Fig. 10 the solution to the Dirac equation in D with zero Dirichlet boundary condition by using the real-space quantum lattice Bolztmann method proposed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF] and illustrating the total wave reflection at the domain boundary when PMLs or ABCs are not used. This simple example illustrates the fact without using PML and or even low-order ABC a wavefunction is totally reflected. This is a standard issue when solving the Dirac equation, in particular when studying laser-molecule interaction. We report in Fig. 11 (left), the maximum of the density as a function of time for σ 0 = 10 -3 , 10 -2 , 10 -1 , 1, 5 in logscale, with the PML of Type IV for θ = π/4. We clearly see that σ 0 = 10 -2 provides the best absorption properties. We next compare in Fig. 12, for σ 0 = 10 -2 , the reflection for θ = π/16, θ = π/8, θ = π/4 and θ = π/2 showing the importance of properly selecting θ, with best absorption at θ = π/16. We report in Fig. 11 (right), the maximum of the density as a function of the time iteration for different values of θ = π/16, π/8, π/4, π/2 with Type IV-PML. We again see an optimal performance with θ = π/16 but θ = π/8 or θ = π/4 also provides some good results, meaning that the PML is relatively stable with θ. In the last test, we compare in Fig. 13, the efficiency of the PML for four different absorbing functions (Types I to IV), with respectively σ(ν) = σ 0 (ν + δ ν ) 2 , σ(ν) = σ 0 (ν + δ ν ) 3 , σ(ν) = -σ 0 /ν and σ(ν) = σ 0 /ν 2 . We see that the best results are obtained for the Type-IV PML, a similar quality being also obtained for the Types V and VI PMLs (not reported here).

These tests clearly show the relevance of the combination of the TSSP method with PMLs, as we simultaneously benefit from the efficiency and accuracy of FFTs and the high absorption feature of PMLs to avoid reflections at the boundary.