Luc Robbiano
email: luc.robbiano@uvsq.fr

Counting function for interior transmission eigenvalues

Keywords: Interior transmission eigenvalues, Weyl law, AMS 2010 subject classification: 35P10, 35P20, 35J57

In this paper we give results on the counting function associated with the interior transmission eigenvalues. For a complex refraction index we estimate of the counting function by Ct n . In the case where the refraction index is positive we give an equivalent of the counting function.

Introduction

In this paper we give an estimate to the counting function associated with the interior transmission eigenvalues. We recall the problem. Let Ω be a smooth bounded domain in R n . Let n(x) be a smooth function defined in Ω, called the refraction index. We say that k = 0 is a interior transmission eigenvalue if there exists (w, v) = (0, 0) such that

       ∆w + k 2 n(x)w = 0 in Ω, ∆v + k 2 v = 0 in Ω, w = v on ∂Ω, ∂ ν w = ∂ ν v on ∂Ω, (1)
where ∂ ν is the exterior normal derivative to ∂Ω. We consider here the function n(x) complex valued. In physical models, we have n(x) = n 1 (x) + in 2 (x)/k where n j are real valued. Taking u = w -v and ṽ = k 2 v, we obtain the following equivalent system if k = 0,

   ∆ + k 2 (1 + m) u + mv = 0 in Ω, (∆ + k 2 )v = 0 in Ω, u = ∂ ν u = 0 on ∂Ω, (2)
where, for simplicity, we have replaced ṽ by v and n by 1 + m.

When k ∈ R, this problem is related with scattering problem. We can find a precise result in Colton and Kress [START_REF] David | Inverse acoustic and electromagnetic scattering theory[END_REF]Theorem 8.9] first proved by Colton, Kirsch and Päivärinta [START_REF] David | Far-field patterns for acoustic waves in an inhomogeneous medium[END_REF] and in a survey by Cakoni and Haddar [START_REF] Fioralba | Transmission Eigenvalues in Inverse Scattering Theory[END_REF].

As the problem si not self-adjoint even for n(x) real valued, usual tools used in self-adjoint cases cannot be applied, in particular, even for operator with compact resolvent, the existence of k is not always true.

A lot of results was obtained this last years using several methods. When n(x) is real, Päivärinta and Sylvester [START_REF] Lassi | Transmission eigenvalues[END_REF] proved that there exist interior transmission eigenvalues; Cakoni, Gintides, and Haddar [START_REF] Fioralba | The existence of an infinite discrete set of transmission eigenvalues[END_REF] proved that the set of k 2 j is infinite and discrete. For n(x) complex valued Sylvester [START_REF] John | Discreteness of Transmission Eigenvalues via Upper Triangular Compact Operators[END_REF] proved that this set is discrete finite or infinite. In [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF] we proved that there exist a infinite number of complex eigenvalues and the associated generalized eigenspaces span a dense space in L 2 (Ω) ⊕ L 2 (Ω).

Other related problems are studied in literature, problem with cavities studied by Cakoni, Çayören, and Colton [START_REF] Fioralba | Transmission eigenvalues and the nondestructive testing of dielectrics[END_REF], Cakoni, Colton, and Haddar [START_REF] Fioralba | The interior transmission problem for regions with cavities[END_REF], problem for operators of order m > 2 by Hitrik, Krupchyk, Ola, and Päivärinta [START_REF] Michael | Transmission eigenvalues for operators with constant coefficients[END_REF][START_REF] Michael | Transmission eigenvalues for elliptic operators[END_REF].

Lakshtanov and Vainberg [START_REF] Evgeny | Ellipticity in the Interior Transmission Problem in Anisotropic Media[END_REF][START_REF] Evgeny | Remarks on interior transmission eigenvalues, Weyl formula and branching billiards[END_REF][START_REF] Evgeny | Bounds on positive interior transmission eigenvalues[END_REF] studied the counting function for problems with different boundary conditions. For Problem [START_REF] Shmuel | Lectures on elliptic boundary value problems[END_REF], in [START_REF] Evgeny | Applications of elliptic operator theory to the isotropic interior transmission eigenvalue problem[END_REF] they obtain a lower estimate of the counting function for real interior transmission eigenvalues.

For counting function in [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF] we gave some non optimal estimate. This estimate was improved by Dimassi and Petkov [START_REF] Mouez | Upper bound for the counting function of interior transmission eigenvalues[END_REF]. Their estimate have the same size than the one found below in Theorem 8 except for a factor 3 √ 3. For constant m there is a recent result given by Pham and Stefanov [START_REF] Ha | Weyl asymptotics of the transmission eigenvalues for a constant index of refraction[END_REF] where they give an equivalent of counting function in this case.

The main results of this paper are Theorem 7 and Theorem 8. In Theorem 7 we prove than the counting function satisfies an estimate in Ct n and we prove that if we denote by λ j the eigenvalues of the problem

j∈N 1 λ p j -z p = (2π) -n |z| -p+n/2 Ω ((1 + m(x) -p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx+o(|z| -p+n/2 ),
when |z| goes to ∞ and z in a line outside a domain related with the range of n(x).

In the case where n(x) > 0 is real valued or if n 1 (x) > 0 in the case n(x) = n 1 (x) + n 2 (x)/k this estimate allows, applying a tauberian theorem to give an equivalent of the counting function. We find that N (t) ∼ αt n where the precise value of α is given in Theorem 8. These results are also proven by Faierman [START_REF] Melvin | Transmission eigenvalues for parameter-elliptic boundary problems[END_REF] in a preprint. The methods used are very close that the one used here but he assumes that n(x) = 1 every where. This condition excludes the case of cavity. Here we assume only that n(x) = 1 in a neighborhood of the boundary.

Notations and background

Let Ω be a C ∞ bounded domain in R n . Let n(x) ∈ C ∞ (Ω) be complex valued. We set m(x) = n(x) -1. We consider also the case where n(x) = n 1 (x) + in 2 (x)/k where n j (x) are real valued and k the spectral parameter. This case is different of the previous one but can be treated similarly. We assume that for all x ∈ Ω, n(x) = 0, or n 1 (x) = 0 or equivalently m(x) = -1. We assume that there exists a neighborhood W of ∂Ω such that for x ∈ W , n(x) = 1 or n 1 (x) = 1 or equivalently m(x) = 0. Actually if n(x) = 1 for all x ∈ ∂Ω, such a neighborhood W exists.

We denote by C e the cone in C defined by

C e = {z ∈ C, ∃x ∈ Ω, ∃λ ≥ 0, such that z = λ(1 + m(x))}. (3)
In the case where n

(x) = n 1 (x) + in 2 (x)/k, C e = [0, ∞) if n 1 (x) > 0 for all x ∈ Ω, and C e = [-∞, 0] if n 1 (x) < 0 for all x ∈ Ω.
Here we give some notations useful for the statement of the results. We use the notations and the results proven in [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF], except some change of sign.

Let z ∈ C, we denote by B z (u, v) = (f, g) the mapping defined from

H 2 0 (Ω) ⊕ {v ∈ L 2 (Ω), ∆v ∈ L 2 (Ω)} to L 2 (Ω) ⊕ L 2 (Ω) by -1 1+m ∆ -z u -m 1+m v = f in Ω (-∆ -z)v = g in Ω (4)
In the case where n = n 1 +in 2 /k we must change the definition of B z . We define m 1 (x) = n 1 (x)-1 and m 2 (x) = n 2 (x). The mapping Bk (u, v) = (f, g) is given by

-1 1+m1 ∆ -k 2 -ik m2 1+m1 u + -m1 1+m1 -im2 k(1+m1) v = f in Ω (-∆ -k 2 )v = g in Ω (5)
Remark that the principal symbol of Bk is the same than the one of B z if we set z = k 2 . Under an assumption on z, B z is invertible for some z.

Theorem 1. Assume C e = C, then there exists z ∈ C such that B z is a bijective map from

H 2 0 (Ω) ⊕ {v ∈ L 2 (Ω), ∆v ∈ L 2 (Ω)} to L 2 (Ω) ⊕ L 2 (Ω).
In the case n(x) = n 1 (x) + n 2 (x)/k, here C e = C and we have the same result.

Theorem 2. There exists k ∈ C such that Bk is bijective from H 2 0 (Ω) ⊕ {v ∈ L 2 (Ω), ∆v ∈ L 2 (Ω)} to L 2 (Ω) ⊕ L 2 (Ω). If for z ∈ C the solution B z (u, v) = (f, g) exists we denote by R z (f, g) = (u, v). In case where n(x) = n 1 (x) + in 2 (x)/k if for k ∈ C the solution of Bk (u, v) = (f, g) exists we denote by Rk (f, g) = (u, v). Theorem 3. Assume C e = C, there exists z ∈ C such that the resolvent R z from H 2 (Ω) ⊕ L 2 (Ω) to itself is compact.
In particular, applying the Riesz theory, the spectrum is finite or is a discrete countable set. If λ = 0 is in the spectrum, λ is an eigenvalue associated with a finite dimensional generalized eigenspace.

Theorem 4. There exists k ∈ C such that the resolvent Rk from H 2 (Ω) ⊕ L 2 (Ω) to itself is compact.

In particular, we can apply the Riesz theory, the spectrum is finite or a discrete countable set. If λ = 0 is in the spectrum, λ is an eigenvalue associated with a finite dimensional generalized eigenspace.

Remark 1. Actually if z 0 ∈ C e ∪[0, ∞) for all λ > 0 large enough we can take z = λz 0 in Theorems 1 and 3.

If k 2 0 ∈ C e ∪ [0, ∞) for all λ > 0 large enough we can take k = λk 0 in the Theorems 2 and 4. Here we estimate the resolvent in the exterior of a conic neighborhood of C e ∪ [0, ∞). In particular if n 1 (x) > 0, the eigenvalues k 2 are in all small conic neighborhood of (0, +∞), except for a finite number of eigenvalues.

In general for a non self-adjoint problem, we cannot claim that the spectrum is non empty. In the following theorem, with a stronger assumption on C e , we can prove that the spectrum is non empty.

We say that C e is contained in a sector with angle less than θ if there exist θ 1 < θ 2 , such that C e ⊂ {z ∈ C, z = 0 or z |z| = e iϕ , where θ 1 ≤ ϕ ≤ θ 2 }, and θ 2 -θ 1 ≤ θ.

Theorem 5. Assume that C e is contained in a sector with angle less than θ with θ < 2π/p where 4p > n and θ < π/2. Then there exists z such that the spectrum of R z is infinite and the space spanned by the generalized eigenspaces is dense in

H 2 0 (Ω) ⊕ {v ∈ L 2 (Ω), ∆v ∈ L 2 (Ω)}.
Theorem 6. There exists k such that the spectrum of Rk is infinite and the space spanned by the generalized eigenspaces is dense in

H 2 0 (Ω) ⊕ {v ∈ L 2 (Ω), ∆v ∈ L 2 (Ω)}.
Remark 2. These results are based on the theory given in Agmon [START_REF] Shmuel | Lectures on elliptic boundary value problems[END_REF] and using the spectral results on Hilbert-Schmidt operators. In this theory we deduce that the spectrum is infinite from the proof that the generalized eigenspaces form a dense subspace in the closure of the range of R z [resp. Rk ].

In [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF] we proved that

R p z [resp. Rp k ] is a Hilbert-Schmidt operator if 4p > n. We can deduce the spectral decomposition of R z [resp. Rk ] from that of R p z [resp. Rp k ].
Let z j be the elements of the spectrum of R z [resp. Rk ] and E j the generalized associated eigenspace. We denote by

N (t) = |zj | -1 ≤t 2 dim E j .
If z j is an eigenvalue of R z , λ j = -z + 1/z j is an eigenvalue of B 0 and we have N (t) ∼ ♯{j, |λ j | ≤ t 2 }, where λ j are counted with multiplicity.

Results

We denote by ω j for j = 1, • • • , p, the roots of z p = 1.

Theorem 7. We assume as in Theorem 5 that θ < 2π/p and θ < π/2 where p satisfies 2p > n and 4p > 4 + n. Then, there exists C > 0 such that N (t) ≤ Ct n .

Moreover let µ ∈ C such that |µ| = 1 and we assume that

ω j µ ∈ C e ∪ (0, +∞) for j = 1, • • • , p. We denote by a(x) = (1 + m(x)) -1 or a(x) = (1 + m 1 (x)) -1 if n(x) = n 1 (x) + k -1 n 2 (x).
We fix z 0 such that the resolvent R z0 exist and let µ j such that 1/µ j are the eigenvalues of R z0 counted with multiplicity. Then we have

j∈N 1 µ p j -z p = (2π) -n |z| -p+n/2 Ω (a p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx + o(|z| -p+n/2 ), (6)
when z = rµ and r goes to ∞.

Remark 3. The first part of the theorem improve [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF]Theorem 7] where we found the estimate

N (t) ≤ Ct n+4 . Theorem 8. We assume n 1 (x) = 1 + m 1 (x) > 0 for all x ∈ Ω when n(x) = n 1 (x) + n 2 /k. Then N (t) = αt n + o(t n ) where α = (2π) -n Vol(B 1 ) Ω ((1 + m 1 (x)) n/2 (x) + 1)dx.
Remark 4. By a more precise study in a neighborhood of the boundary we can obtain a result with a smaller remainder in (6) but this estimate does not allow to prove a result better on the counting function. Malliavin [START_REF] Malliavin | Un théorème taubérien relié aux estimations de valeurs propres[END_REF] was proved a tauberian theorem with a sharp remainder but this requires an estimate on j∈N 1 λ p j -z p in a complex domain except in a parabola neighborhood of (0, ∞). Here the estimate are proved in a complex domain except in a conic neighborhood of (0, ∞). It is maybe possible to improve this result following Hitrik, Krupchyk, Ola and Päivärinta [START_REF] Michael | The interior transmission problem and bounds on transmission eigenvalues[END_REF] where they prove that the eigenvalues are in a parabolic neighborhood of (0, ∞).

Proof of Theorem 7

The proof is based on Lemmas 4.1 and 4.2 below. We introduce some notations.

We set S = R z0 and T = S p where p satisfies the assumption of Theorem 7. We set

T λ = T (I -λT ) -1 .
As T λ is a matrix of operators we denote

T λ = T 11 T 12 T 21 T 22 ,
We denote by

K 11 K 12 K 21 K 22 , (7)
the kernel of T z p .

Lemma 4.1. Under the assumption of Theorem 7, there exists C > 0 such that N (t) ≤ Ct n . Moreover let V a conical neighborhood of C e ∪ [0, ∞), there exists R > 0 such all z ∈ C, satisfying |z| ≥ R and ω j z / ∈ V , we have

Ω K 11 (x, x)dx + Ω K 22 (x, x)dx = j∈N 1 λ p j -z p . ( 8
)
Lemma 4.2. With the notation of Theorem 7 we have

|z| p-n/2 Ω K 11 (x, x)dx + Ω K 22 (x, x)dx goes to (2π) -n Ω (a p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx when |z| goes to ∞.
Clearly Lemmas 4.1 and 4.2 imply Theorem 7.

Proof of Lemma 4.1

We recall that ω j for j = 1, • • • , p, are the roots of z p = 1, we have

(1 -z p S p ) = p j=1 (1 -ω j zS). (9)
The operator (1 -z p S p ) is invertible if and only if (1 -ω j zS) is invertible for all j. Thus we have

T z p = S p p j=1 (1 -ω j zS) -1 = p j=1 S ωj z .
If we denote by

S z = S 11 S 12 S 21 S 22 .
We recall that S = R z0 and We apply the results of [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF]. The estimates below are given by [22, theorem 10] for k ≥ 1. The estimate on S 12 and S 22 for k = 0 are also given by [22, theorem 10]. For k = 0, [22, Proposition 2.3] gives the estimate on S 21 and [22, Lemma 2.1] gives the estimate on S 11 . The relation between z and h is zh 2 = µ, in particular h 2 |z| = 1.

S z = (R z0 ) z = R z0 (I -zR z0 ) -1 = (R -1 z0 -z) -1 = (B z0 -z) -1 = (B 0 -z 0 -z) -1 = R z0+z , (10)
S 11 : H 2k sc (Ω) → H 2k+2 sc (Ω) with S 11 H 2k sc (Ω)→H 2k+2 sc (Ω) ≤ C|z| -1 S 12 : H 2k sc (Ω) → H 2k+4 sc (Ω) with S 12 H 2k sc (Ω)→H 2k+4 sc (Ω) ≤ C|z| -2 S 21 : H 2k sc (Ω) → H 2k sc (Ω) with S 21 H 2k sc (Ω)→H 2k sc (Ω) ≤ C S 22 : H 2k sc (Ω) → H 2k+2 sc (Ω) with S 22 H 2k sc (Ω)→H 2k+2 sc (Ω) ≤ C|z| -1 (11)
We denote by

Λ z = |z| 0 0 1/ |z|.
We remark that

Λ z AΛ -1 z = A 11 |z|A 12 (1/|z|)A 21 A 22 , where A = A 11 A 12 A 21 A 22 . (12)
We deduce from [START_REF] Michael | Transmission eigenvalues for operators with constant coefficients[END_REF] that

Λ z S ωj z Λ -1 z : L 2 (Ω) ⊕ L 2 (Ω) → H 2 sc (Ω) ⊕ L 2 (Ω) with an operator norm less than C|z| -1 and Λ z S ωjz Λ -1 z : H 2k+2 sc (Ω) ⊕ H 2k sc (Ω) → H 2k+4 sc (Ω) ⊕ H 2k+2 sc
(Ω) with an operator norm less than C|z| -1 .

As

Λ z T z p Λ -1 z = p j=1 Λ z S ωj z Λ -1 z , (13)
we deduce that

Λ z T z p Λ -1 z : L 2 (Ω) ⊕ L 2 (Ω) → H 2p sc (Ω) ⊕ H 2p-2 sc (Ω), (14)
with an operator norm less than C|z| -p . We can prove that N (t) ≤ Ct n . First we weaken [START_REF] Karamata | Neuer Beweis und Verallgemeinerung einiger Tauberian-Sätze[END_REF] to consider

Λ z T z p Λ -1 z as a map between L 2 (Ω) ⊕ L 2 (Ω) → H 2p-2 sc (Ω) ⊕ H 2p-2 sc
(Ω), with an operator norm less than

C|z| -p . As v H s ≤ h -2s v H s sc , we obtain Λ z T z p Λ -1 z L 2 (Ω)⊕L 2 (Ω)→H 2p-2 (Ω)⊕H 2p-2 (Ω) ≤ C|z| -1 , Λ z T z p Λ -1 z L 2 (Ω)⊕L 2 (Ω)→L 2 (Ω)⊕L 2 (Ω) ≤ C|z| -p . (15)
We can apply the theorem 13.5 in Agmon [START_REF] Shmuel | Lectures on elliptic boundary value problems[END_REF], that is, if m > n/2 we have

|||T ||| ≤ C T n/(2m) m T 1-n/(2m) 0
, where |||T ||| is the Hilbert-Schmidt norm and T m is the operator norm of the map

L 2 (Ω)⊕L 2 (Ω) → H m (Ω) ⊕ H m (Ω).
We apply this estimate with m = 2p -2 > n/2 and we have

|||Λ z T z p Λ -1 z ||| |z| -n 4(p-1) |z| -(1-n 4(p-1) )p = |z| -p+n/4 .
We can follow the proof of Theorem 7 in [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF]. If we denote by µ j complex numbers such that µ -1 j are eigenvalue of S counted with multiplicity, then 1 µ p j -z p are the eigenvalues of T z p and thus the eigenvalues of Λ z T z p Λ -1 z . We obtain

j 1 |µ p j -z p | 2 ≤ |||Λ z T z p Λ -1 z ||| 2 ≤ C|z| -2p+n/2 . Let µ ∈ C such that |µ| = 1, and ω j µ ∈ C e ∪ (0, ∞) for all j = 1, • • • , p. We take z = t 2 µ. If |µ j | ≤ t 2 , we have |µ p j -z p | ≤ 2t 2p . Then we have |µj |≤t 2 1 4t 4p ≤ j 1 |µ p j -z p | 2 ≤ |||T z p ||| 2 ≤ Ct -4p+n .
Then we obtain N (t) ≤ Ct n . Now we prove Formula (6). Estimate [START_REF] Karamata | Neuer Beweis und Verallgemeinerung einiger Tauberian-Sätze[END_REF] implies that

T 11 L 2 (Ω)→H 2p sc (Ω) ≤ C|z| -p and as v H s ≤ h -2s v H s sc , we have T 11 L 2 (Ω)→H 2p (Ω) ≤ C. ( 16
)
To estimate the norm of T 22 , we shall use that S 12 is a mapping from L 2 (Ω) to H 4 sc (Ω). Actually if we take g ∈ L 2 (Ω), we have

Λ z S ωj z Λ -1 z (0, g) ∈ H 4 sc (Ω)⊕ H 2 sc (Ω) and Λ z S ωj z Λ -1 z (0, g) H 4 sc (Ω)⊕H 2 sc (Ω) ≤ C|z| -1 g L 2 (Ω)
. We can repeat the previous argument for the p -1 other factors Λ z S ωj z Λ -1 z and we obtain that

Λ z T z p Λ -1 z (0, g) H 2p+2 sc (Ω)⊕H 2p sc (Ω) ≤ C|z| -p g L 2 (Ω) .
In particular this means that T 22 g H 2p sc (Ω) ≤ C|z| -p g L 2 (Ω) , and

T 22 g H 2p (Ω) ≤ C g L 2 (Ω) . (17)
By [START_REF] Evgeny | Ellipticity in the Interior Transmission Problem in Anisotropic Media[END_REF], we have

T 11 L 2 (Ω)→L 2 (Ω) + T 22 L 2 (Ω)→L 2 (Ω) ≤ C|z| -p . ( 18
)
We can apply the theorem 13.9, Agmon [START_REF] Shmuel | Lectures on elliptic boundary value problems[END_REF]. If 2p > n we have for j = 1, 2, K jj ∈ H q (Ω × Ω)

where q = 2p -[n/2] -1 > n/2.
In particular the trace K jj (x, x) is well defined in L 2 (Ω) and

|K jj (x, x)| 2 dx 1/2 ≤ C T jj L 2 (Ω)→H 2p (Ω) + T * jj L 2 (Ω)→H 2p (Ω) n/(2p) T jj 1-n/(2p) L 2 (Ω)→L 2 (Ω) + T jj L 2 (Ω)→L 2 (Ω) (19)
Remark 5. The adjoint of B z is given by an analogous formula than (4). Indeed we find that the adjoint B * z (p, q) = (g 1 , g 2 ) is given by

-∆((1 + m) -1 p) -zp = g 1 in Ω -∆q -zq + m(1 + m) -1 p = g 2 in Ω q |∂Ω = ∂ ν q |∂Ω = 0 on ∂Ω.
Using the relation between R z and S z (see ( 10)), we deduce that the adjoint of S z satisfies the same estimate than S z given in [START_REF] Michael | Transmission eigenvalues for operators with constant coefficients[END_REF]. By [START_REF] Michael | The interior transmission problem and bounds on transmission eigenvalues[END_REF], the adjoint of T z satisfies the same properties than T z given in ( 16), ( 17) and [START_REF] Evgeny | Applications of elliptic operator theory to the isotropic interior transmission eigenvalue problem[END_REF].

For j = 1, 2,(19) implies, from ( 16), ( 17) and ( 18)

Ω |K jj (x, x)| 2 dx 1/2 ≤ C|z| -p+(n/2) (20)
We recall that λ -1 j are the eigenvalues of S counted with multiplicity the eigenvalues of T are λ -p j . The indices are such that |λ j | ≤ |λ j+1 |. As N (t) ≤ Ct n , this implies that |λ j | ≥ Cj 2/n where C > 0.

In particular 1/|λ j | p converges if 2p > n. By Theorem 12.17 in Agmon [START_REF] Shmuel | Lectures on elliptic boundary value problems[END_REF], there exists a constant c ∈ C such that

Tr(T T z p ) = j∈N 1 (µ p j -z p )µ p j + c. (21)
We recall that the trace is defined in the theorem 12.20 in Agmon [START_REF] Shmuel | Lectures on elliptic boundary value problems[END_REF] for an operator

Q = Q 1 Q 2 where Q 1 and Q 2 are Hilbert-Schmidt operators. Moreover if K is the kernel of Q, K(x, x) is definite for almost all x, we have Ω |K(x, x)|dx < ∞ and
Tr(Q) = Ω K(x, x)dx. (22)
We remark as by assumption µ p j and z p are not in the same cone, we have

|µ p j -z p | ∼ |µ p j | + |z p |. Then 1 |µ p j -z p ||µ p j | ≤ C |µ p j | and 1 |µ p j -z p ||µ p j | → 0 when |z| → +∞.
This implies that j∈N 1 (µ p j -z p )µ p j → 0 when |z| → +∞.

In [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF] we proved that |||T z p ||| ≤ C|z| 1-p+n/4 as 1 -p + n/4 < 0 thus |||T z p ||| goes to 0 as |z| goes to +∞. We have | Tr(T T z p )| ≤ |||T ||||||T z p ||| goes to 0 as |z| goes to +∞. Then c = 0 in [START_REF] Ha | Weyl asymptotics of the transmission eigenvalues for a constant index of refraction[END_REF]. We obtain that

Tr(z p T T z p ) = j∈N 1 µ p j -z p - 1 µ p j .
We have z p T T z p = z p T 2 (I -z p T ) -1 = -T (I -z p T )(I -z p T ) -1 + T (I -z p T ) -1 = T z p -T . By Formula ( 22) the trace of z p T T z p is given by the integral of its kernel and, as the integral of the kernel of T z p exists by [START_REF] Lassi | Transmission eigenvalues[END_REF], the integral of trace of kernel of z p T T z p -T z p does not depend of z, we obtain

Ω K 11 (x, x)dx + Ω K 22 (x, x)dx = j∈N 1 µ p j -z p - 1 µ p j + c (23)
By [START_REF] Lassi | Transmission eigenvalues[END_REF], for j = 1, 2, Ω K jj (x, x)dx goes to 0 when |z| goes to 0 and as

1 |µ p j -z p | ≤ C 1

Proof of Lemma 4.2

We recall some facts on pseudo-differential operators. Let a(x, ξ) be in C ∞ (R n × R n ) we say that a is a symbol of order m if for all α, β ∈ N n , there exist C α,β > 0, such that

|∂ α x ∂ β ξ a(x, ξ)| ≤ C α,β ξ m-|β| ,
where ξ 2 = 1 + |ξ| 2 . In particular a polynomial in ξ of order m with coefficients in C ∞ (R n ) with bounded derivatives of all orders, is a symbol of order m.

With a symbol we can associate an semi-classical operator by the following formula

Op(a)u = a(x, D)u = 1 (2π) n e ixξ a(x, hξ)û(ξ)dξ = 1 (2hπ) n e ixξ/h a(x, ξ)û(ξ/h)dξ.
If a(x, ξ) is a symbol of order m (which can depend of h), if a(x, ξ) = b(x, ξ) + hc(x, ξ) where b(x, ξ) and c(x, ξ) are symbols of order m, we call b(x, ξ) the principal symbol of a(x, ξ) which is definite modulo h. This formula makes sense for u ∈ S (R n ) and we can extend it to u ∈ H s for all s. For a, a symbol of order m, there exists C > 0 such that for all u ∈ H s , a(x, D)u H s-m sc ≤ C u H s sc . In the following, when we use pseudo-differential operator we have always cut-off functions supported in Ω in each side of the operator. We do not have to consider the action of pseudo-differential operator on H s (Ω) space as in [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF].

We begin with a description on S z in all compact set in Ω.

Lemma 4.3. Let µ ∈ C such that |µ| = 1 and we assume that µ ∈ C e ∪ (0, +∞)

for j = 1, • • • , p. Let z ∈ C such that z/µ ∈ (0, ∞) is large enough. Let θ and θ be functions in C ∞ 0 (Ω) such that θ(x) = 1 if x in the support of θ, and θ(x) = 1 if x in a compact subset of Ω. Then we have θΛ z S z Λ -1 z = |z| -1 θBθ + |z| -1/2 W θΛ z S z Λ -1 z , ( 24
)
where W = Λ z KΛ -1 z and W * are bounded on H s (Ω) and the principal symbol of B is given on the support of θ by

(a|ξ| 2 -µ) -1 (a|ξ| 2 -µ) -1 V (|ξ| 2 -µ) -1 0 (|ξ| 2 -µ) -1 .
Proof. We apply the result proved in [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF]. Let us recall the notations and the main results. We multiply Equations ( 4) by h 2 , we denote by µ = h 2 z where µ belongs to a bounded domain of C, a = 1/(1 + m) and V = m/(1 + m). We change (f, g) in (-f, -g).

We recall the assumption made on m, we have m(x) = -1 for all x ∈ Ω and m(x) = 0 for x in a neighborhood of ∂Ω.

Thus following (4), we obtain the system

   -ah 2 ∆ -µ u -h 2 V v = h 2 f in Ω (-h 2 ∆ -µ)v = h 2 g in Ω u = ∂ ν u = 0 on ∂Ω. (25)
In case where n(x) = n 1 (x) + n 2 (x)/k we multiply Equations ( 5) by h 2 , we denote by µ = -h 2 k 2 , a = 1/(1 + m 1 ), V = m 1 /(1 + m 1 ) + hm 2 /(ν + νm 1 ) where ν = hk. We change (f, g) in (-f, -g). Thus following ( 5) we obtain the system

   -ah 2 ∆ + hW 0 -µ u -h 2 V v = h 2 f in Ω (-h 2 ∆ -µ)v = h 2 g in Ω u = ∂ ν u = 0 on ∂Ω, (26)
where

W 0 = -νm 2 /(1 + m 2 ). In particular the principal semi-classical symbol of -ah 2 ∆ + hW 0 -µ is -a|ξ| 2 -µ, the principal semi-classical symbol of V is m 1 /(1 + m 1 )
. In what follows only the principal symbols of -ah 2 ∆ + hW 0 -µ and -h 2 ∆ -µ must be take account. For simplicity we write the proof for the system (25), the case of system (26) may be treat following the same way. Now we compute the symbol of the resolvent in Ω. Let φ 0 , φ 1 and φ 2 ∈ C ∞ 0 (Ω) where φ 0 φ 1 = φ 0 and φ 1 φ 2 = φ 1 . We take φ 0 such that φ 0 = 1 on the support of θ. Let Q be a parametrix of -ah 2 ∆ -µ such that φ 0 Qφ 1 (-ah 2 ∆ -µ) = φ 0 -hK where K is of order -1 and the principal symbol of Q is (a|ξ| 2 -µ) -1 . As

φ 1 -ah 2 ∆ -µ φ 2 u = h 2 V φ 1 v + h 2 φ 1 f in R n , applying φ 0 Q to this equation, we obtain φ 0 u = h 2 φ 0 Qφ 1 f + h 2 φ 0 Q (V φ 1 v) + hKφ 2 u. (27)
We apply the same method use above to compute v (see [START_REF] Luc | Spectral analysis on interior transmission eigenvalues[END_REF]Section 2.2]). Let φ 3 ∈ C ∞ 0 (Ω) where φ 2 φ 3 = φ 2 and θ = 1 on the support of φ 3 . The choices of the φ j are compatible with θ and θ. We have φ 1 Qφ 2 (-h 2 ∆ -µ) = φ 1 -hK -1 where K -1 is of order -1 and the principal symbol of

Q is (|ξ| 2 -µ) -1 .
We apply φ 2 on the equation on v in (25), we have

φ 2 (-h 2 ∆ -µ)φ 3 v = h 2 φ 2 g.
Applying the parametrix φ 1 Q, we have

φ 1 v = h 2 φ 1 Qφ 2 g + hK -1 φ 3 v. (28)
With this equation and ( 27) we obtain

φ 0 u = h 2 φ 0 Qφ 1 f + h 4 φ 0 Q V φ 1 Qφ 2 g + hKφ 2 u + h 3 K-1 φ 3 v, (29)
where K-1 is an operator of order -1 . Let

A = h 2 Qφ 1 h 4 QV φ 1 Qφ 2 0 h 2 φ 1 Qφ 2 and K = K 11 h 2 K 12 0 K 22 ,
where

K 11 = Kφ 2 , K 12 = K-1 and K 22 = K -1 φ 2 .
In particular K jk is bounded on H s sc (Ω) for all s ≥ 0. Indeed, all the operators contain cut-off thus

K jk φ 3 u is compactly supported in Ω if u ∈ H s sc (Ω).
We recall that S = R z0 and by [START_REF] Melvin | Transmission eigenvalues for parameter-elliptic boundary problems[END_REF], S z = R z0+z , if the resolvent R z0+z exists. In what follows, z 0 is fixed and we have the relation z = -µ/h 2 -z 0 . In particular |z| -1/2 ∼ h for large |z|. With these relations we have, following (28) and (29), S z (f, g) = (u, v), and

φ 0 S z = φ 0 Aφ 3 + |z| -1/2 Kφ 3 S z .
Thus we can write

φ 0 Λ z S z Λ -1 z = |z| -1 φ 0 Bφ 3 + |z| -1/2 W φ 3 Λ z S z Λ -1 z , (30)
where W = Λ z KΛ -1 z is bounded on H s (Ω) and the principal symbol of B is given on the support of φ 0 by

(a|ξ| 2 -µ) -1 (a|ξ| 2 -µ) -1 V (|ξ| 2 -µ) -1 0 (|ξ| 2 -µ) -1 .
As θφ 0 = θ and φ 3 θ = φ 3 , (30) gives (24). As K is a semi-classical pseudo-differential operator, W * is also bounded on H s (Ω) Remark 6. Formula (24) does not give a description on the operator S z in Ω. It gives only S z in all compact in Ω. In the proof below we need also estimates on S z up the boundary given in [START_REF] Michael | Transmission eigenvalues for operators with constant coefficients[END_REF] to absorb the error terms.

Lemma 4.3 gives the principal symbol of S z , the following lemma gives the principal symbol of T z p . Lemma 4.4. Let p ∈ N \ {0}. Let ϕ 0 and ϕ 1 be functions in C ∞ 0 (Ω) such that ϕ 0 = 1 on a neighborhood of supp ϕ 1 .

ϕ 1 Λ z T z p Λ -1 z ϕ 0 = |z| -p ϕ 1 B p ϕ 0 + |z| -p-1/2 ϕ 1 R p ϕ 0 , (31)
where ϕ 1 R p ϕ 0 satisfies the following property if we denote by

ϕ 1 R p ϕ 0 = R 11 p R 12 p R 21 p R 22 p R 11 p : L 2 (Ω) → H 2p sc (Ω) R 12 p : L 2 (Ω) → H 2p+2 sc (Ω) R 21 p : L 2 (Ω) → H 2p-2 sc (Ω) R 22 p : L 2 (Ω) → H 2p sc (Ω), (32)
where the norm operator are uniformly bounded with respect h, and the principal symbol of B p is

(a p |ξ| 2p -µ p ) -1 Q -2p-2 (x, ξ) 0 (|ξ| 2p -µ p ) -1 (33)
Moreover the adjoint of R jq p satisfies (32) where the norm operator are uniformly bounded with respect h.

Proof. We argue by induction on k and for that we must introduce a sequence of cut-off functions. Let χ k and χk be cut-off functions such that χk χk+1 = χk+1 , χ k χk = χk , χk χ k+1 = χ k+1 . We can assume that χ0 = 1 on the support of ϕ 0 and and χp = 1 on the support of ϕ 1 . We can apply Formula (24) where θ is replaced by χk and θ by χ k . We have

χk Λ z S z Λ -1 z = |z| -1 χk Bχ k + |z| -1/2 W χ k Λ z S z Λ -1 z , (34)
where W and W * are bounded on H s (Ω).

Now we can finish the proof of Lemma 4.2. We take ϕ 0 such that ϕ 0

(x) = 1 is d(x, R n \ Ω) ≥ 2δ and ϕ 0 (x) = 0 is d(x, R n \ Ω) ≤ δ. We take ϕ 1 such that ϕ 1 (x) = 1 is d(x, R n \ Ω) ≥ 4δ and ϕ 1 (x) = 0 is d(x, R n \ Ω) ≤ 3δ.
With the notation of (31) we deduce from (32) that for j = 1, 2 we have

|z| -p-1/2 ϕ 1 R jj p ϕ 0 L 2 (Ω)→H 2p sc (Ω) |z| -p-1/2 .
This implies as for T jj in ( 17) and ( 18) that

|z| -p-1/2 ϕ 1 R jj p ϕ 0 L 2 (Ω)→L 2 (Ω) |z| -p-1/2 and |z| -p-1/2 ϕ 1 R jj p ϕ 0 L 2 (Ω)→H 2p (Ω) |z| -1/2 .
By Formula (19) applied to the kernel of |z| -p-1/2 ϕ 1 R jj p ϕ 0 denoted by K R jj and the properties on R p and its adjoint given in Lemma 4.3, we obtain

Ω |K R jj (x, x)| 2 dx 1/2 ≤ C|z| -p+(n/2)-1/2 . ( 38
)
By the principal symbol of B p given in Lemma 4.4 we can compute the integral of the kernel K jj (x, x). Denoting by b(x, ξ) either (a p |ξ| 2p -µ p ) -1 or (|ξ| 2p -µ p ) -1 , the kernel of a diagonal term is given by (2πh) -n e i(x-y)ξ/h ϕ 1 (x)b(x, ξ)ϕ 0 (y)dξ and this integral make sense if p > n/2. Denoting by K Bp jj (x, y) the diagonal terms of the kernel of |z| -p ϕ 1 B p ϕ 0 , we obtain

Ω K Bp 11 (x, x)dx + Ω K Bp 22 (x, x)dx = (2π) -n |z| -p+n/2 ϕ 1 (x) (a p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx + O(|z| -p+(n+1)/2 ), (39)
where the error term O(|z| (n+1)/2 ) is given by the lower order terms in the symbolic calculus. From (31), ( 38) and (39), we deduce

ϕ 0 (x)(K 11 (x, x) + K 22 (x, x))dx = (2π) -n |z| -p+n/2 ϕ 1 (x) (a p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx + O(|z| -p+(n+1)/2
).

Now we can write

Λ z T z p Λ -1 z = ϕ 1 Λ z T z p Λ -1 z ϕ 0 + ϕ 1 Λ z T z p Λ -1 z (1 -ϕ 0 ) + (1 -ϕ 1 )Λ z T z p Λ -1 z . If K(x, y) is the kernel of Λ z T z p Λ -1
z , the kernel of the left hand side terms are respectively, ϕ 1 (x)K(x, y)ϕ 0 (y), ϕ 1 (x)K(x, y)(1 -ϕ 0 (y)) and (1 -ϕ 1 (x))K(x, y).

In particular, from the properties of the supports of ϕ j , we have ϕ

1 (x)K(x, x)(1 -ϕ 0 (x)) = 0. Let F (x) = K 11 (x, x) + K 22 (x, x), to show that Ω |z| p-n/2 F (x)dx -(2π) -n Ω (a p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx ≤ Cδ + C δ | z| -1/2 , (40) We shall prove |z| p-n/2 ϕ 1 (x)F (x)dx -(2π) -n Ω (a p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx ≤ Cδ + C δ | z| -1/2 , (41) and |z| p-n/2 Ω (1 -ϕ 1 (x))F (x)dx ≤ Cδ 1/2 . (42)
Obviously ( 41) and ( 42) imply (40), and (40) implies Lemma 4.2. To prove (41), we apply (39). We have

{x, d(x,R n \Ω)≤2δ} (a p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx ≤ C {x, d(x,R n \Ω)≤2δ} ξ -2p dξdx ≤ Cδ,
indeed |ξ| 2p and a p |ξ| 2p are not in the same cone as µ p by assumption. Thus the principal term in (41) given by the principal term from (39), can be estimate by Cδ and the error terms from (39) can be estimate by C δ |z| -1/2 . To prove (42), using [START_REF] Lassi | Transmission eigenvalues[END_REF] we obtain

Ω (1 -ϕ 1 (x))K jj (x)dx ≤ Ω (K jj (x, x)) 2 dx 1/2 {x,d(x,R n \Ω)≤4δ} dx 1/2 ≤ C|z| -p+n/2 δ 1/2 ,
which is Estimate (42).

Proof of Theorem 8

We shall apply the following tauberian Theorem cited in Agmon [1, th. 14.5], The proof is given in Karamata [START_REF] Karamata | Neuer Beweis und Verallgemeinerung einiger Tauberian-Sätze[END_REF].

Theorem 9 (Tauberian Theorem). Let σ(λ) be a non decreasing function for λ > 0, let 0 < a < 1, let α be a non-negative number, and suppose that as t → +∞, +∞ 0 dσ(λ) λ + t = αt a-1 + o(t a-1 ).

Then as λ → +∞ σ(λ) = α sin πa πa λ a + o(λ a ).

Theorem 8 is implied by tauberian theorem and Theorem 7 which gives

j∈N 1 µ p j -z p = A|z| -p+n/2 + o(|z| -p+n/2 ), where A = (2π) -n Ω (a p |ξ| 2p -µ p ) -1 + (|ξ| 2p -µ p ) -1 dξdx, with µ = z/|z|. (43)
We take z = t 1/p e iπ/p , where t > 0, we obtain j∈N 1 µ p j + t = At -1+n/2p + o(t -1+n/2p ).

Let µ j = δ j + iν j where δ j and ν j are real. Let A(t) = j∈N 1 δ p j +t . By assumption (see Remark 1), for all ε > 0, there exists R ε > 0, such that if δ j ≥ R ε then |ν j | ≤ εδ j .

We have for j large enough such that |ν j | ≤ εδ j ,

1 δ p j + t - 1 µ p j + t ≤ Cε 1 δ p j + t , (45)
where C depends only on p.

Denoting N ε such that for all j ≥ N ε , then |ν j | ≤ εδ j . We have

We recall that µ = z/|z| and z = t 1/p e -iπ/p thus µ = -1. To compute A from (43), we must compute integral as (b p |ξ| 2p + 1) -1 dξ, where b > 0. We have

For all ε > 0, there exists J ε such that for all j ≥ J ε , |ν j | ≤ εδ j then |µ j | ≤ (1 + ε)δ j . Thus there exist C ε > 0 such that ♯{j ∈ N, δ j ≤ t 2 } ≤ ♯{j ∈ N, (1 + ε) -1 |µ j | ≤ t 2 } + C ε , which is equivalent to ♯{j ∈ N, δ j ≤ (1 + ε) -1 t 2 } ≤ ♯{j ∈ N, |µ j | ≤ t 2 } + C ε .

We obtain

(1 + ε) -n/2 α ≤ lim inf t→∞ t -n ♯{j ∈ N, |µ j | ≤ t 2 }.
As this estimate is true for all ε > 0 and from (49) we have

α ≤ lim inf t→∞ t -n ♯{j ∈ N, |µ j | ≤ t 2 } ≤ lim sup t→∞ t -n ♯{j ∈ N, |µ j | ≤ t 2 } ≤ α.
This is the statement of Theorem 8.

  p | → 0 when |z| → +∞, we obtain c = j∈N 1 µ p j . Then (23) gives the statement of Lemma 4.1.

  (b p |ξ| 2p + 1) -1 dξ = nVol(B 1 ) +∞ 0 r n-1 (b p r 2p + 1) -1 dr = n(2p) -1 b -n/2 Vol(B 1 ) +∞ 0 σ n/(2p)-1 (σ + 1) -1 dσ = n(2p) -1 b -n/2 Vol(B 1 )π sin -1 (nπ/(2p)),where the last integral is computed by residue theorem (see Cartan[6, p 107]).This gives that2pA sin nπ/(2p) nπ = (2π) -n Vol(B 1)Ω (a -n/2 (x) + 1)dx = α.Now to prove the statement of Theorem 8, we must prove that α = lim t -n ♯{j ∈ N, δ j ≤ t 2 } = lim t -n ♯{j ∈ N, |µ j | ≤ t 2 } when t goes to ∞. Except for a finite number of values, δ j > 0 and asδ j ≤ |µ j | there exists C > 0 such that ♯{j ∈ N, |µ j | ≤ t 2 } ≤ ♯{j ∈ N, δ j ≤ t 2 } + C.Thus lim sup t→∞ t -n ♯{j ∈ N, |µ j | ≤ t 2 } ≤ α.

  + O(N ε /t) + O(εA(t)) = At -1+n/2p + o(t -1+n/2p ). (47)For ε small enough we deduce there exists C > 0 such that(1/C)t -1+n/2p ≤ A(t) ≤ Ct -1+n/2p . At -1+n/2p + o(t -1+n/2p ).Denoting σ(λ) = ♯{j ∈ N, δ p j ≤ λ} where the number is counted with multiplicity. We have At -1+n/2p + o(t -1+n/2p ).

	Nε j=1	1 j + t| |µ p	≤ C	Nε j=1	1 j + t δ p	≤ CN ε /t.	(46)
	We deduce from (44), (45) and (46) that
	A(t) Using this in (47) we obtain
	j∈N = j∈N 1 δ p j + t 1 δ p j + t = +∞ 0 dσ(λ) λ + t = By tauberian Theorem we obtain
	♯{j ∈ N, δ p j ≤ λ} = 2pA	sin nπ/(2p) nπ	λ

n/(2p) + o(λ n/(2p) ), which is equivalent to ♯{j ∈ N, δ j ≤ t 2 } = 2pA sin nπ/(2p) nπ t n + o(t n ).

We prove by recurrence the following formula

where the semi-classical principal symbol of B k is given by

where Q -2k-2 (x, ξ) is a symbol of order -2k -2. The operators R jq k and their adjoints satisfy Estimates (32) with p = k.

For k = 1, Formula (36) and properties (32) for R jq 1 and their adjoints follow from ( 12), [START_REF] Michael | Transmission eigenvalues for operators with constant coefficients[END_REF], Remark 5 and (34).

If Formula (35) is true for k we have

where

By (34) for k + 1 and (35) for k, we have

The term |z| -1-k χk+1 Bχ k+1 χk B k χ0 gives the first right hand side term of (35), where B k+1 = Bχ k+1 χk B k and the principal symbol is given by Formula (36) on the support of χk+1 .

The three other terms have the form of R k+1 and satisfy Estimates (32). Indeed, the power of |z| is obtained as the operator norm of Λ z S z Λ -1 z is bounded by |z| -1 . To prove the mapping between the H s sc , we denote by A q a generic operator of order q mapping H s to H s-q . We check that

The properties on adjoints follow from the recurrence assumptions on R k , the properties on W * and Remark 5. By (34) for k + 1 and as χ k+1 (1 -χk ) = 0, we have

By [START_REF] Michael | Transmission eigenvalues for elliptic operators[END_REF] and [START_REF] Michael | Transmission eigenvalues for operators with constant coefficients[END_REF] the operator norm of this term is |z| -k-1/2 . The proof that L 2 satisfies Estimates (32) for p = k + 1 is obtained by (37). The properties on the adjoint of L 2 follow from Remark 5 and the properties on W * . From (36) for k = p, and as z p -µ p = p j=1 (z -ω j µ), we obtain that p j=1 (a|ξ| 2 -ω j µ) -1 = (a p |ξ| 2p -µ p ) -1 . This gives the diagonal terms of the symbol of B p in Formula (33).


Sébastien Boyaval
email: sebastien.boyaval@enpc.fr

Marco Picasso
email: marco.picasso@epfl.ch

A posteriori analysis of Chorin-Temam scheme for Stokes equations

Keywords: Operator splitting, Pressure corrrection, Projection method, A posteriori error estimation French: Séparation d'opérateurs, Correction de pression, Méthode

We consider Chorin-Temam scheme (the simplest pressure-correction projection method) for the time-discretization of an unstationary Stokes problem in D ⊂ R d (d = 2, 3) given µ, f : (P) find (u, p) solution to u| t=0 = u 0 , u| ∂D = 0 and

Inspired by the analyses of the Backward Euler scheme performed by C.Bernardi and R.Verfürth, we derive a posteriori estimators for the error on ∇u in L 2 (0, T ; L 2 (D))-norm. Our invesigation is supported by numerical experiments.

French version: On discrétise en temps par le schéma Chorin-Temam un problème de Stokes non-stationnaire posé dans D ⊂ R d (d = 2, 3) étant donnés µ, f : (P) trouver (u, p) solution de u| t=0 = u 0 , u| ∂D = 0 et (1). En s'insipirant des analyses de C.Bernardi and R.Verfürth pour le schéma Euler rétrograde, nous construisons des estimateurs a posteriori pour l'erreur commise sur ∇u en norme L 2 (0, T ; L 2 (D)). Notre étude est étayée par des expériences numériques.

donc l'estimateur (22) plutôt que (20) (tiré directement de Prop. 2 et Prop. 3) si on veut une estimation robuste (c'està-dire de qualité indépendante des paramètres de discrétisation). D'autre part, bien que notre estimation ne soit pas totalement efficace (comme dans [START_REF] Kharrat | Time error estimators for the chorin-temam scheme[END_REF][START_REF] Kharrat | A posteriori error analysis of time-dependent stokes problem by chorin-temam scheme[END_REF], nos estimateurs ne sont pas bornés inférieurement et supérieurement par l'erreur), on montre néanmoins numériquement qu'elle peut être utile dans certains cas, et en particulier qu'elle est plus précise que celle proposée dans [START_REF] Kharrat | Time error estimators for the chorin-temam scheme[END_REF][START_REF] Kharrat | A posteriori error analysis of time-dependent stokes problem by chorin-temam scheme[END_REF] (plus de termes d'erreur sont pris en compte).

Pour des approximations (à λ > 0 donné) des composantes du vecteur vitesse et de la pression u = π sin(λt) sin(2πy) sin(πx) 2 ;sin(2πx) sin(πy) 2 p = sin(λt) cos(πx) sin(πy) avec des éléments finis continus P 2 et P 1 par morceaux dans D ≡ (-1, 1) × (-1, 1) (d = 2) maillé régulièrement avec des simplexes, on a calculé numériquement l'efficacité des estimateurs (20), ( 22) et (23) pour t ∈ (0, T ) discrétisé avec des pas de temps constants ∆t = T/N (N ∈ N). En effet, notre analyse a posteriori du cas semi-discret en temps se prolonge au cas complètement discret (en décomposant les résidus discrets en composantes temporelles et spatiales comme dans [START_REF] Verfürth | A posteriori error analysis of space-time finite element discretizations of the time-dependent Stokes equations[END_REF] on obtient directement les versions discrètes en espace des estimateurs semi-discrets en temps plus des estimateurs pour l'erreur en espace), et l'erreur de discrétisation en espace est par ailleurs négligeable ici pour notre exemple numérique (comme observé dans [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] où il est utilisé pour λ = 1). 23)) donc on ne le voit que pour ∆t assez petit même si on ajoute le terme div u ∆t 2 L ∞ (0,t;Q d ) à l'estimateur (22). Sans parler de l'estimation de l'erreur sur u en norme L ∞ (0, T ; L 2 (D)), on n'a donc pas encore totalement résolu le problème de trouver un estimateur efficace et robuste pour l'erreur commise sur ∇u en norme L 2 (0, T ; L 2 (D)) par le schéma Chorin-Temam. Il faudrait au moins ajouter des coefficients devant les termes de l'estimateur (22) plus div u ∆t 2 L ∞ (0,t;Q d ) si on veut l'utiliser en pratique. Néanmoins, nous espérons que cette étude apporte un nouvel éclairage à la question. the usual L 2 inner-products for scalar and vector functions in D and introduce the standard functional spaces [START_REF] Temam | Navier-Stokes equations[END_REF][START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF]

Q := L 2 (D) , Q := {q ∈ L 2 (D) , D q = 0} , W := [H 1 0 (D)] d , V := {v ∈ [H 1 0 (D)] d , div v = 0} . (2)
We consider a weak formulation of problem (P) with µ > 0, f ∈ L 2 (0, T ; Q d ) (given in a Bochner space), u 0 ∈ V: (P) find u ∈ L 2 (0, T ; W) and p ∈ L 2 (0, T ; Q) such that u(0) = u 0 in V, and the following equation holds in L 2 (0, T )

d dt (u, v) + µ(∇u, ∇v) -(p, div v) + (q, div u) = ( f , v) , ∀(v, q) ∈ W × Q . (3)
It is well-known that problem (P) is well-posed [START_REF] Temam | Navier-Stokes equations[END_REF][START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF] (in particular, u ∈ C([0, T ], V) so initial condition makes sense) and because of the regularity assumptions, it also holds

∂ t u ∈ L 2 ((0, T ) × D), p ∈ L 2 (0, T ; H 1 (D)) and in L 2 (0, T ) (∂ t u, v) + µ(∇u, ∇v) + (∇p, v) -(∇q, u) = ( f , v) , ∀(v, q) ∈ W × H 1 (D) . (4)
A standard time-discretization of (4) is Chorin-Temam scheme [START_REF] Chorin | Numerical Solution of the Navier-Stokes Equations[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF]: given u -1/2 = u 0 , p 0 = 0, for n = 0 . . . N -1,

given ∆t n ∈ (0, ∆t], f ∆t = 1 ∆t n t n+1 t n f (s)ds (t n = n-1 k=0 ∆t k ; t N = T ), (P n ) find u n+1/2 ∈ W, p n+1 ∈ Q ∩ H 1 (D) solutions to ( u n+1/2 -u n-1/2 ∆t n + ∇p n , v) + µ(∇u n+1/2 , ∇v) = ( f n+1 , v) ∀v ∈ V , (5a)
1 ∆t n+1 (div u n+1/2 , q) = -(∇p n+1 , ∇q) ∀q ∈ Q , (5b)
which yields approximations whose rate of convergence to solutions of (P) is well-known a priori [START_REF] Guermond | Some implementations of projection methods for Navier-Stokes equations[END_REF][START_REF] Prohl | Projection and quasi-compressibility methods for solving the incompressible Navier-Stokes equations[END_REF][START_REF] Guermond | On stability and convergence of projection methods based on pressure Poisson equation[END_REF]:

Proposition 1.
The following estimate holds:

u ∆t -u L 2 (0,T ;W) + p ∆t -p L 2 (0,T ;Q) = O(∆t 1 2 ) as ∆t → 0 , (6)
where u ∆t and p ∆t are defined as

u ∆t (t) = t -t n ∆t n u n+1/2 - t -t n+1 ∆t n u n-1/2 , p ∆t (t) = p n . ∀t ∈ (t n , t n+1 ] , (7)
In this work, we would like to numerically evaluate a posteriori the time discretization error with a view to adequately choosing the time steps ∆t n of Chorin-Temam scheme in practice (under a given error tolerance), which is still an open problem. A posteriori error estimations have been proposed for the Backward-Euler scheme (including full discretizations, in time and space) [1, 2] but they do not straightforwardly apply here. And a posteriori analyses of Chorin-Temam scheme have indeed been carried out recently [START_REF] Kharrat | Time error estimators for the chorin-temam scheme[END_REF][START_REF] Kharrat | A posteriori error analysis of time-dependent stokes problem by chorin-temam scheme[END_REF] but they suggest an estimator (different than ours) that does not account for the whole error. The present invesigation focuses on fully computable error bounds for Chorin-Temam scheme derived from the generic a posteriori framework introduced in [START_REF] Verfürth | A posteriori error analysis of space-time finite element discretizations of the time-dependent Stokes equations[END_REF] for the unstationary Stokes equations. Although our estimator is a priori not fully efficient, it is better than other ones and useful in some cases.

Note that in the following, we denote by a b any relation a ≤ Cb between two real numbers a, b where C > 0 is a numerical constant independent of the data of the problem. Moreover, we shall use standard inequalities such as

div v Q ≤ d 1/2 ∇v Q d×d , ∀v ∈ W and Poincaré-Friedrichs inequality with constant C P (D) > 0, then also max( v 2 Q d , ∇v 2 Q d×d ) ≤ v 2 W ≤ (1 + C 2 P ) ∇v 2 Q d×d , ∀v ∈ W . (8)
In Section 2, we derive a posteriori error estimates following the procedure of [START_REF] Verfürth | A posteriori error analysis of space-time finite element discretizations of the time-dependent Stokes equations[END_REF], i.e. invoking Π : W → W, a projection such that v -Πv ∈ V. For all v ∈ W, Πv is the solution of Stokes equations:

∃! q v ∈ Q, ∃Υ(D) > 0 such that (∇Πv, ∇w) = (q v , div w) (r, div Πv) = (r, div v) , ∀(w, r) ∈ W × Q , (9a)
Υ ∇Πv Q d×d ≤ div v Q . ( 9b
)
In Section 3, we numerically test our a posteriori estimator.

A posteriori estimation of semi-discrete errors

Let us define residuals for Chorin-Temam approximations u ∆t , p ∆t as in [START_REF] Guermond | Some implementations of projection methods for Navier-Stokes equations[END_REF] of the solution to the problem (P)

< R u , v > W ,W = ( f , v) -(∂ t u ∆t , v) -(∇p ∆t , v) -µ(∇u ∆t , ∇v) ≡ ( f -f ∆t , v) + µ(∇u ∆t,+ -∇u ∆t , ∇v) , ∀v ∈ W , (10a) (R p , q) = -(div u ∆t , q) , ∀q ∈ Q , ( 10b
)
where

f ∆t = f n+1 , u ∆t,+ = u n+ 1 2 for t ∈ (t n , t n+1 ].
The errors e u = uu ∆t , e p = pp ∆t satisfy:

(∂ t e u + ∇e p , v) + µ(∇e u , ∇v) + (div e u , q) =< R u , v > W ,W +(R p , q) , ∀(v, q) ∈ W × Q . ( 11
)
Testing [START_REF] Kharrat | A posteriori error analysis of time-dependent stokes problem by chorin-temam scheme[END_REF] against v = e u -Πe u ∈ V, q = 0, yields in D (0, T ) (distributional sense)

1 2 d dt e u 2 Q d + µ ∇e u 2 Q d×d =< R u , e u > W ,W -< R u , Πe u > W ,W +(∂ t e u , Πe u ) + µ(∇e u , ∇Πe u ) . ( 12
)
Using Young inequality with [START_REF] Prohl | Projection and quasi-compressibility methods for solving the incompressible Navier-Stokes equations[END_REF], div e u = R p , (9b), Πe u = -Πu ∆t , e u (0) = 0, and integrating by part, one obtains

e u 2 L ∞ (0,t;Q d ) + µ ∇e u 2 L 2 (0,t;Q d×d ) R u 2 L 2 (0,t;W ) + µ R p 2 L 2 (0,t;Q) + t 0 |(e u , ∂ t Πe u )| + (e u , Πe u ) L ∞ (0,t) . (13)
If we follow [START_REF] Verfürth | A posteriori error analysis of space-time finite element discretizations of the time-dependent Stokes equations[END_REF], then (13) yields a computable upper-bound using the following inequalities with Young's one

t 0 (e u , ∂ t Πu ∆t ) ≤ e u L ∞ (0,t;Q d ) Π∂ t u ∆t L 1 (0,t;Q d ) (14a) (e u , Πu ∆t ) L ∞ (0,t) e u L ∞ (0,t;Q d ) div u ∆t L ∞ (0,t;Q d ) . (14b) Since 1 ∂ t e u + ∇e p 2 L 2 (0,T ;W ) ≤ 2 R u 2 L 2 (0,T ;W ) + 2 ∇e u 2
L 2 (0,T ;Q d×d ) also holds from (10a), one indeed obtains from (9b):

Proposition 2. There exists a constant c + (D) > 0 such that the following computable estimations hold

1 c + max e u 2 L ∞ (0,T ;Q d ) , ∂ t e u + ∇e p 2 L 2 (0,T ;W ) , µ ∇e u 2 L 2 (0,T ;Q d×d ) ≤ f -f ∆t 2 L 2 (0,T ;Q d ) + µ ∇u ∆t,+ -∇u ∆t 2 L 2 (0,T ;Q d×d ) + µ div u ∆t 2 L 2 (0,T ;Q) + div ∂ t u ∆t 2 L 1 (0,T ;Q) + div u ∆t 2 L ∞ (0,T ;Q) . (15)
On the other hand, from (10a), (10b) and div u = 0, one has

µ ∇u ∆t,+ -∇u ∆t 2 L 2 (0,T ;Q d×d ) + µ div u ∆t 2 L 2 (0,T ;Q) f -f ∆t 2 L 2 (0,T ;Q d ) + ∂ t e u + ∇e p 2 L 2 (0,T ;W ) + µ ∇e u 2 L 2 (0,T ;Q d×d ) , (16)
from which one next straightforwardly obtains the counterpart of (15) if one uses, in addition to (16),

div ∂ t u ∆t 2 L 1 (0,T ;Q) T min n=0...N-1 |∆t n | 2 ∇e u 2 L 2 (0,T ;Q d×d ) . (17)
Proposition 3. There exists a constant c -(D) > 0 such that the following computable lower bound holds

c -µ ∇u ∆t,+ -∇u ∆t 2 L 2 (0,T ;Q d×d ) + µ div u ∆t 2 L 2 (0,T ;Q) + 1 N div ∂ t u ∆t 2 L 1 (0,T ;Q) ≤ f -f ∆t 2 L 2 (0,T ;Q d ) + ∂ t e u + ∇e p 2 L 2 (0,T ;W ) + µ + 1 min n=0...N-1 |∆t n | ∇e u 2 L 2 (0,T ;Q d×d ) . ( 18
)
Proof of (17). We use the following inequality with div u = 0, noting 6

∆t n t n+1 t n div u ∆t 2 Q ≥ div u n+ 1 2 2 Q + div u n-1 2 2 Q : div ∂ t u ∆t 2 L 1 (0,T ;Q) ≤ N N-1 n=0 div(u n+ 1 2 -u n-1 2 ) 2 Q ≤ 2N N-1 n=0 ( div u n+ 1 2 2 Q + div u n-1 2 2 Q ) ≤ N-1 n=0 12N ∆t n t n+1 t n div u ∆t 2 Q .
So the framework introduced in [START_REF] Verfürth | A posteriori error analysis of space-time finite element discretizations of the time-dependent Stokes equations[END_REF] for an a posteriori analysis of a Backward Euler discretization of Stokes problem still applies here with Chorin-Temam scheme (it applies with any scheme provided the reconstructions u ∆t , p ∆t are defined using appropriate discrete variables). Though, the point is now to let not only the residuals, but also the two last terms in (13), be easily and sharply estimated (contrary to the fully discrete Backward Euler case in [START_REF] Verfürth | A posteriori error analysis of space-time finite element discretizations of the time-dependent Stokes equations[END_REF], these terms cannot be neglected here because they can be of the same order as the error). We draw the following conclusions. First, Prop. 2 and 3 suggest that the procedure of [START_REF] Verfürth | A posteriori error analysis of space-time finite element discretizations of the time-dependent Stokes equations[END_REF] should be modified here to estimate the error

µ ∇e u 2 L 2 (0,T ;Q d×d ) + ∂ t e u + ∇e p 2 L 2 (0,T ;W ) (19)
a posteriori in a more robust way than by the estimator (20) obtained straightforwardly from the estimations above:

µ ∇u ∆t,+ -∇u ∆t 2 L 2 (0,T ;Q d×d ) + µ div u ∆t 2 L 2 (0,t;Q) + div ∂ t u ∆t 2 L 1 (0,T ;Q) . ( 20
)
1 Observe that the convergence of ∂ t u ∆t + ∇p ∆t to ∂ t u + ∇p in L 2 (0, T ; W ) is natural here, like for Backward-Euler schemes [START_REF] Bernardi | A posteriori error analysis of the fully discretized time-dependent Stokes equations[END_REF].

4

For instance, if one replaces (14a) with the following upper bound (21), on noting ( 8) and (9b),

t 0 |(e u , ∂ t Πe u )| ∇e u L 2 (0,t;Q d ) div ∂ t u ∆t L 2 (0,t;Q) , (21)
then bounds similar to (15) and (18) hold but with div ∂ t u ∆t L 2 (0,t;Q) instead of div ∂ t u ∆t L 1 (0,t;Q) and without invoking discretization parameters like N and ∆t n , which suggests the a posteriori error estimator (22) more robust than (20):

µ ∇u ∆t,+ -∇u ∆t 2 L 2 (0,T ;Q d×d ) + µ div u ∆t 2 L 2 (0,T ;Q) + div ∂ t u ∆t 2 L 2 (0,T ;Q) . (22)
Of course, this is not a fully efficient estimator yet, since it is a priori not bounded above and below by the error (19), even if one neglects the source error ff ∆t2 L 2 (0,T ;Q d ) of "high" order O(∆t 2 ) -recall ( 6) -. It is nevertheless useful in some cases, as shown in the next section. Second, (22) sometimes improves some estimators in the literature like

µ ∇u ∆t,+ -∇u ∆t 2 L 2 (0,T ;Q d×d ) + N-1 n=0 ∆t n+1 ∇p n+1 -∆t n ∇p n 2 Q (23)
that was proposed in [START_REF] Kharrat | Time error estimators for the chorin-temam scheme[END_REF][START_REF] Kharrat | A posteriori error analysis of time-dependent stokes problem by chorin-temam scheme[END_REF]. Clearly, for small ∆t, our estimator is larger than the one proposed in [START_REF] Kharrat | Time error estimators for the chorin-temam scheme[END_REF][START_REF] Kharrat | A posteriori error analysis of time-dependent stokes problem by chorin-temam scheme[END_REF], on noting

∆t n+1 ∇p n+1 -∆t n ∇p n 2 Q div u n+1/2 -div u n-1/2 2 Q ∆t ∆t n div ∂ t u ∆t 2 Q (t)
for t ∈ (t n , t n+1 ], after using a Poincaré inequality with (5b). And, the numerical example of the following Section 3 indeed shows that ( 22) is a better upper-bound than (23), at least when the error is not mainly driven by div u ∆t 2 L ∞ (0,T ;Q) .

Numerical results

We want to bring numerical evidences that estimator ( 22) is sometimes i) useful and ii) better than (20) and ( 23). Given λ > 0, we numerically compute the efficiencies of the three estimators using discrete approximations of u = π sin(λt) sin(2πy) sin(πx) 2 ;sin(2πx) sin(πy) 2 p = sin(λt) cos(πx) sin(πy) in D ≡ (-1, 1) × (-1, 1) (d = 2), with t ∈ (0, T ) uniformly discretized by time steps ∆t = T/N (N ∈ N) when µ = 1. We discretize in space the velocity components and the pressure with, respectively, continuous P 2 and P 1 Finite-Elements functions, i.e. in conforming discrete spaces W h ⊂ W, Q h ⊂ (Q ∩ H 1 (D)) defined on regular simplicial meshes of D. The a posteriori analysis of Section 2 still applies with right-hand side in (13) defined using now fully-discrete approximations. Then indeed, following [START_REF] Verfürth | A posteriori error analysis of space-time finite element discretizations of the time-dependent Stokes equations[END_REF], one can decompose the fully-discrete residuals in a sum of two terms, one accounting for space-discretization errors and one for time-discretization errors. The two last terms in the (new) right-hand side of (13) remain the same (they are explicitly computable). This yields estimators linked to the time discretization which are exactly the space-discrete counterparts of the terms in the bounds (15) and (18). Moreover, in our numerical example, space discretization errors prove negligible in comparison with time discretization errors (as already observed in [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for λ = 1). We thus next show only numerical results obtained for one sufficiently fine mesh (with more than 10 5 vertices).

We compare the effectivities of (space-discrete versions of) the a posteriori error estimators (20), ( 22) and (23) evoked in the previous section for the (space-discrete) error ∇e u h 2 L 2 (0,T ;Q d×d ) + ∂ t e u h + ∇e p h 2 L 2 (0,T ;W ) . One clearly sees from the numerical results obtained for λ = 10, T ≤ 3 in Fig. 1 that i) (20) is not robust when ∆t is too small or T too large compared with (22), and ii) ( 22) is better than the estimator (23) in so far as, for that specific case, it has the same decay rate than the error (19) and not a faster one like (23). Though, our estimator (22) is still not fully efficient, even when adding the term div u ∆t 2 L ∞ (0,T ;Q) to (22), since it is not bounded above and below by the error. Furthermore, if we use it as such (that is as a sum of terms without coefficients), in some cases, it also fails (like (23)) at evaluating correctly the error. For instance when λ = 1, the error (19) scales like div u ∆t 2 L ∞ (0,T ;Q) with respect to ∆t, while the other terms in (22) are of higher-order in ∆t. But this cannot be observed unless ∆t is very small, even if we use (22) plus div u ∆t 2 L ∞ (0,T ;Q) as an estimator 2 insofar as the magnitude of the latter term is much smaller than the former (10 -1 vs. 10 2 ). Then, for too large ∆t, the effectivity of our estimator also decays, and the error still cannot be evaluated confidently. So, without even mentionning the error e u L ∞ (0,T ;Q d ) , the question how to estimate a posteriori error discretizations in Chorin-Temam scheme efficiently and robustly (in all cases) thus remains open. One should at least coefficient adequately the terms in the estimator above (22) plus div u ∆t 2 L ∞ (0,T ;Q) . We nevertheless hope to have shed new light on the problem.

1 .

 1 Numerical solutions to Stokes equations by Chorin-Temam pressure-correction projection method Given a smooth bounded open set D ⊂ R d (d = 2, 3) with boundary ∂D of class C 2 , let us denote similarly by (•, •)

Figure 1 :

 1 Figure 1: For ∆t = .1, .05, .025, .0125, .00625, effectivities in log scale (as a function of T ) of (20) -top left -, (22) -top right -, and (23) -bottom left -( div u ∆t 2 L ∞ (0,T ;Q) included) at estimating (19) when λ = 10, T ≤ 3; and div u ∆t 2 L ∞ (0,T ;Q) / error (19) -bottom right -when λ = 1, T ≤ 10.

Note that in fact we also added the term div u ∆t 2 L ∞ (0,T ;Q) to (20), (22) and (23) in Fig.1, but it is small compared to other terms, thus unseen.


Jan Bulla
email: bulla@math.unicaen.fr

Christophe Chesneau
email: chesneau@math.unicaen.fr

Maher Kachour
email: kachour.maher@gmail.com

On the bivariate Skellam distribution

Keywords: Skellam distribution, Bivariate Skellam distribution, Maximum likelihood estimate. 2000 Mathematics Subject Classification 62G30, 62E10

In this paper, we introduce a new distribution on Z 2 , which can be viewed as a natural bivariate extension of the Skellam distribution. The main feature of this distribution a possible dependence of the univariate components, both following univariate Skellam distributions. We explore various properties of the distribution and investigate the estimation of the unknown parameters via the method of moments and maximum likelihood. In the experimental section, we illustrate our theory. First, we compare the performance of the estimators by means of a simulation study. In the second part, we present two applications to a real data set and show how an improved fit can be achieved by estimating mixture distributions.

Introduction

In recent years, there has been a growing interest in studying bivariate discrete random variables, in an effort to explain phenomena in various areas of application. More specifically, bivariate discrete random variables are appropriate for modelling paired count data exhibiting correlation. Paired count data occur in many contexts such as marketing (joint purchases of two products), econometrics (number of voluntary and involuntary job changes), insurance (number of accidents in a site before and after infrastructure changes), engineering (faults due to different causes), epidemiology (joint concurrence of two different diseases), and many others.

The Bivariate Poisson Distribution (BPD), originally derived by [START_REF] Mckendric | Applications of mathematics to medical problems[END_REF] as a solution to a differential equation arising in a biological application, is probably the best known bivariate discrete distribution. For a comprehensive treatment of the bivariate Poisson distribution and its multivariate extensions, the reader can refer to the books of Kocherlakota and Kocherlakota (1992) and [START_REF] Johnson | Multivariate Discrete Distributions[END_REF], and the review articles by [START_REF] Papageorgiou | Multivariate discrete distributions[END_REF] and [START_REF] Kocherlakota | Bivariate discrete distributions[END_REF]. The BPD has found many applications recently (e.g., [START_REF] Karlis | Analysis of Sports Data by Using Bivariate Poisson Models[END_REF] have used it for modelling sports data). Many other bivariate discrete distributions have been introduced, for details we refer to [START_REF] Marshall | Multivarirate distribution generated from mixtures of convolution and product families[END_REF], [START_REF] Boucher | Models of insurance claim counts with time dependence based on generalization of poisson and negative binomial distributions[END_REF], and [START_REF] Cheon | Tests for independence in a bivariate negative binomial model[END_REF]. In particular, we cite the Bivariate Negative Binomial Distribution (BNBD), which was used in insurance theory for describing of the number of accidents in transportation, and the Bivariate Geometric Distribution (BGD), which is considered plausible for reliability modelling. Note that the literature on bivariate count distributions with negative correlation is limited. One of the reasons is that negative correlation in bivariate counts occurs rather infrequently. However there are such models in the literature, as for example the bivariate Poisson-lognormal model of [START_REF] Aitchinson | The multivariate Poisson-log normal distribution[END_REF], the finite mixture model developed in [START_REF] Karlis | Finite multivariate Poisson mixtures with applications[END_REF] and models based on copulas see, e.g. [START_REF] Nikoloulopoulos | Finite normal mixture copulas for multivariate discrete data modeling[END_REF] and the references therein.

In contrast to the well-known situation when the paired data are counts, i.e., observed on N 2 , sometimes the data take values in Z 2 . For example, when analysing intra-daily stock prices the changes take both positive and negative integer values, known also as ticks (the price can go up or down on certain predefined ranges of value). The price change is therefore characterized by discrete jumps. While bivariate discrete distribution for non-negative paired data are now abundant, to our best knowledge there is a shortage of bivariate discrete distribution defined on Z 2 . We intend to contribute to this literature by presenting a bivariate discrete distribution based on the Skellam distribution.

The remainder of the paper proceeds as follows. Firstly, we present the considered bivariate Skellam distribution and investigate some of its properties in Section 2. Section 3 describes two possible approaches for estimating the unknown parameters efficiently. Section 4 presents results from a simulation study and a real data set, respectively. Last, Section 5 is devoted to the proofs.

2 The BSkellam(λ 0 , λ 1 , λ 2 ) distribution
Before analyzing the bivariate setting, we recall the definition of the univariate Skellam distribution. Skellam(λ 1 , λ 2 ), if and only if

X d = U 1 -U 2 ,
where U 1 and U 2 are two independent random variables such that U i ∼ Poisson(λ i ) for any i ∈ {1, 2}. The probability mass function of X is given by

P(X = x) = e -(λ1+λ2) λ x 1 ∞ i=max(0,-x) (λ 1 λ 2 ) i (x + i)!i! , x ∈ Z.
The distribution of the difference between two independent Poisson random variables was derived by [START_REF] Irwin | The frequency distribution of the difference between two independent variates following the same Poisson distribution[END_REF] for the case of equal parameters. [START_REF] Skellam | The frequency distribution of the difference between two poisson variates belonging to different populations[END_REF] and [START_REF] Prékopa | On composed Poisson distributions. IV. Remarks on the theory of differential processes[END_REF] discussed the case of unequal parameters. Further details can be found in e.g. [START_REF] Karlis | Bayesian modeling of football outcomes: Using the skellam distribution for the goal difference[END_REF] and [START_REF] Al-Zaid | On the poisson difference distribution inference and applications[END_REF].

Remark 1 The definition of the Skellam distribution can be extended to more than the simple difference of two independent Poisson distributions. Indeed, let X 1 and X 2 be two independent Poisson random variables with parameters θ 1 and θ 2 respectively. Let

Y i = X i + W , for i = 1, 2, where W is a random variable independent of X 1 and X 2 . Thus, Z = Y 1 -Y 2 = X 1 -X 2 also follows a Skellam(θ 1 , θ 2 ) distribution.
Inspired by the trivariate reduction method derived by [START_REF] Johnson | Distributions in statistics: Discrete distribution[END_REF], the purpose of this paper is to study a natural bivariate extension of the Skellam(λ 1 , λ 2 ) distribution, including dependence between the components. It is described below.

Definition 2 (BSkellam(λ 0 , λ 1 , λ 2 ) distribution) Let λ 0 ≥ 0, λ 1 > 0 and λ 2 > 0. We say that the bivariate random variable (X 1 , X 2 ) follows a bivariate Skellam distribution, denoted by BSkellam(λ 0 , λ 1 , λ 2 ), if and only if

X 1 ∼ Skellam(λ 1 , λ 0 ), and X 2 ∼ Skellam(λ 2 , λ 0 ).
The probability mass function of (X 1 , X 2 ) is given by

P(X 1 = x 1 , X 2 = x 2 ) = e -(λ1+λ2+λ0) λ x1 1 λ x2 2 ∞ i=max(0,-x1,-x2) (λ 0 λ 1 λ 2 ) i (x 1 + i)!(x 2 + i)!i! for all (x 1 , x 2 ) ∈ Z 2 . Remark 2 Note that Z = X 1 -X 2 ∼ Skellam(λ 1 , λ 2 ).
Hence, the probability distribution of Z is independent of λ 0 . Moreover, for λ 0 = 0, the bivariate Skellam distribution reduces to the product of two independent Poisson distributions.

Lemma 1 below investigates some of the basic properties of the BSkellam(λ 0 , λ 1 , λ 2 ) distribution.

Lemma 1 Let λ 0 ≥ 0, λ 1 > 0, λ 2 > 0 and (X 1 , X 2 ) ∼ BSkellam(λ 0 , λ 1 , λ 2 ). Then -the mean of (X 1 , X 2 ) is (λ 1 -λ 0 , λ 2 -λ 0 ) and the covariance matrix equals Σ =       λ 1 + λ 0 λ 0 λ 0 λ 2 + λ 0       . -the correlation coefficient of (X 1 , X 2 ) is cor(X 1 , X 2 ) = λ 0 (λ 1 + λ 0 )(λ 2 + λ 0 ) .
Once again, one can see that λ 0 > 0 implies that X 1 and X 2 are dependent random variables. the characteristic function of (X 1 , X 2 ) can be calculated as φ (X1,X2) (t 1 , t 2 ) = e λ1(e it 1 -1) e λ2(e it 2 -1) e λ0(e -it 1 -it 2 -1) ∀ (t 1 , t 2 ) ∈ Z 2 .

Remark 3 Using the characteristic function of (X 1 , X 2 ) ∼ BSkellam(λ 0 , λ 1 , λ 2 ), the probability mass function of (X 1 , X 2 ) can be expressed by

P (X 1 = x 1 , X 2 = x 2 ) = 1 (2π) 2 π -π π -π φ (X1,X2) (t 1 , t 2 )e -ix1t1 e -ix2t2 dt 1 dt 2 ∀ (x 1 , x 2 ) ∈ Z 2 . (1)
In particular, this yields the following equality

π -π π -π e λ1(e it 1 -1) e λ2(e it 2 -1) e λ0(e -it 1 -it 2 -1) e -ixt1 e -iyt2 dt 1 dt 2 = (2π) 2 e -(λ1+λ2+λ0) λ x 1 λ y 2 ∞ i=0 (λ 0 λ 1 λ 2 ) i (x + i)!(y + i)!i! ∀ (x 1 , x 2 ) ∈ Z 2 . ( 2
)
Lemma 2 below presents a Gaussian approximation result for the BSkellam(λ 0 , λ 1 , λ 2 ) distribution.

Lemma 2 Let λ 0 ≥ 0, λ 1 > λ 0 and λ 2 > λ 0 and suppose that λ 1 and λ 2 are large enough. Then, the distribution of (X 1 , X 2 ) can be approximate by the following bivarirate gaussian distribution

N 2             λ 1 -λ 0 λ 2 -λ 0       ,       λ 1 -λ 0 λ 0 λ 0 λ 2 -λ 0             .
Remark 4 Recall that the distribution of X 1 and X 2 can be approximated by a Gaussian distribution for large values of λ 1 + λ 0 and λ 2 + λ 0 , respectively.

Parameter estimation

This section is devoted to the estimation of the parameters of a BSkellam(λ 0 , λ 1 , λ 2 ) distribution via the method of moments and the method of the maximum likelihood, respectively. Let λ 0 ≥ 0, λ 1 > 0, and λ 2 > 0 be unknown parameters and (X 1,1 , X 2,1 ), . . . , (X 1,n , X 2,n ) be n i.i.d. bivariate random variables with the common distribution BSkellam(λ 0 , λ 1 , λ 2 ).

The method of moments

Set, for any j ∈ {1, 2},

X j = 1 n n i=1 X j,i , and C 1,2 = 1 n -1 n i=1 (X 1,i -X 1 )(X 2,i -X 2 ).
In view of Lemma 1, the method of moment provides the estimator λ0 = C 1,2 for λ 0 , the estimator λj = X j + λ0 for λ j , j ∈ {1, 2}.

Remark 5 Thanks to the strong law of large numbers, note that λj converge a.s. to λ j for j ∈ {0, 1, 2}.

Remark 6 The moment estimates of λ 0 do not exist if C 1,2 < 0. To solve this problem, a modification is done such that the negative estimate of λ 0 is set to zero and the other estimate is set to equal the absolute value of the respective sample mean.

Maximum likelihood estimation

For any z = (x, y) ∈ Z 2 and λ = (λ

0 , λ 1 , λ 2 ) ∈ [0, ∞[ 3 , set G(λ, z) = λ x 1 λ y 2 ∞ i=max(0,-x,-y) (λ 0 λ 1 λ 2 ) i (x + i)!(y + i)!i! . Lemma 3 For any (x, y) ∈ Z 2 holds ∂G(λ, z) ∂λ 1 = x λ 1 G(λ, z) + λ 0 λ 1 G(λ, (x + 1, y + 1)), ∂G(λ, z) ∂λ 2 = y λ 2 G(λ, z) + λ 0 λ 2 G(λ, (x + 1, y + 1)), and
∂G(λ, z) ∂λ 0 = G(λ, (x + 1, y + 1)).
For any i ∈ {1, . . . , n}, set

z i = (x 1,i , x 2,i ) ∈ Z 2 , z = (z 1 , . . . , z n ), and λ = (λ 0 , λ 1 , λ 2 ) ∈ [0, ∞[ 3 . Then, the likelihood function is given by L(λ, z) = n i=1 P(X 1 = x 1,i , X 2 = x 2,i ) = e -n(λ0+λ1+λ2) n i=1 G(λ, z i ).
Therefore, the j-th partial derivative of the log-likelihood function equals

∂ ln L(λ, z) ∂λ j = -n + n i=1 ∂G(λ,zi) ∂λj G(λ, z i ) .
for any j ∈ {0, 1, 2}. From Lemma 3 follows

∂ ln L(λ, z) ∂λ 1 = -n + 1 λ 1 n i=1 x 1,i + λ 0 λ 1 n i=1 G(λ, z i + 1) G(λ, z i ) , ∂ ln L(λ, z) ∂λ 2 = -n + 1 λ 2 n i=1 x 2,i + λ 0 λ 2 n i=1 G(λ, z i + 1) G(λ, z i ) ,
and

∂ ln L(λ, z) ∂λ 0 = -n + n i=1 G(λ, z i + 1) G(λ, z i ) . Let Z i = (X 1,i , X 2,i ) for any i ∈ {1, . . . , n}. The maximum likelihood estimate λ = ( λ0 , λ1 , λ2 ) satisfies                        0 = -n + 1 λ1 n i=1 X 1,i + λ0 λ1 n i=1 G( λ, Z i + 1) G( λ, Z i ) , 0 = -n + 1 λ2 n i=1 X 2,i + λ0 λ2 n i=1 G( λ, Z i + 1) G( λ, Z i ) , 0 = -n + n i=1 G( λ, Z i + 1) G( λ, Z i ) . (3) Therefore, n i=1 G( λ, Z i + 1)/G( λ, Z i ) = n and, a fortiori, λ1 = X 1 + λ0 , λ2 = X 2 + λ0 . (4)
Substituting (4) into the last equation of the system of equations (3), we obtain the following equation for λ0 :

n = n i=1 G((X 1 + λ0 , X 2 + λ0 , λ0 ), Z i + 1) G((X 1 + λ0 , X 2 + λ0 , λ0 ), Z i ) .
Then, λ1 and λ2 follow directly from λ0 and (4), and the maximum likelihood estimates λ are complete. Nevertheless, for practical purposes a direct numerical maximization of expression 1 and 2, respectively, is also a straightforward approach, and applied in the what follows.

Remark 7 The maximum likelihood estimator λ is asymptotically normally distributed N 3 (λ, I -1 (λ)), where I(λ) is the Fisher information matrix with entries

I i,j (λ) = -E ∂ ln L(λ, (X 1 , X 2 )) ∂λ i-1 ∂λ j-1 with (i, j) ∈ {1, 2, 3}.
Traditionally, this matrix can be estimated by I( λ), which allows the constructions of confidence intervals or hypothesis tests.

Extensions of the bivariate Skellam

Having only three parameters, the bivariate Skellam distribution is a rather parsimonious distribution for bivariate data. However, along with the parsimony naturally come certain drawbacks. For example, mean and skewness of the marginal distributions have identical sign, and the standard deviation depends on the same parameters as the mean. In order to increase the flexibility of the bivariate Skellam, we propose two extensions: Firstly, the inclusion of a shift parameter and, secondly, mixtures of bivariate Skellam distributions.

Shifted BSkellam distribution. The inclusion of a shift parameter k = (k 1 , k 2 ) in the bivariate Skellam is straightforward. Then, we obtain a new bivariate random variable (X 1 * , X 2 * ) = (X 1 + k 1 , X 2 + k 2 ) with probability mass P(X 1 * = x 1 , X 2 * = x 2 ) = P(X 1 = x 1 -k 1 , X 2 = x 2 -k 2 ) where (X 1 , X 2 ) ∼ BSkellam(λ 0 , λ 1 , λ 2 ). Basic properties of this distribution are (E(X 1 * ), E(X 2 * )) = (E(X 1 ) + k 1 , E(X 2 * ) + k 2 ) = (λ 1 -λ 0 + k 1 , λ 2 -λ 0 + k 2 )
, and the covariance matrix identical to that of (X 1 , X 2 ). Mixture of BSkellam distributions. Distributions based on the mixture of several bivariate Skellam distributions can be characterized by the probability mass

P(X 1• = x 1 , X 2• = x 2 ) = M m=1 p m P(X 1m = x 1 , X 2m = x 2 ),
where

M m=1 p m = 1, p m ≥ 0 and (X 1m , X 2m ) ∼ BSkellam(λ 0m , λ 1m , λ 2m
) for all m ∈ {1, . . . , M }. Let us observe that, for all i ∈ Z, we have

P(X 1• = i) = p m P(X 1m = i) and P(X 2• = i) = p m P(X 2m = i). It follows that (E(X 1• ), E(X 2• )) = ( M m=1 p m (λ 1m -λ 0m ), M m=1 p m (λ 2m -λ 0m )), cov(X 1• , X 2• ) = M m=1 p m (λ 0m + (λ 1m -λ 0m )(λ 2m -λ 0m )) - M m=1 p m (λ 1m -λ 0m ) M m=1 p m (λ 2m -λ 0m ), V(X 1• ) = M m=1 p m (λ 1m + λ 0m + (λ 1m -λ 0m ) 2 ) - M m=1 p m (λ 1m -λ 0m ) 2 , and
V(X 2• ) = M m=1 p m (λ 2m + λ 0m + (λ 2m -λ 0m ) 2 ) - M m=1 p m (λ 2m -λ 0m ) 2 .
Moments of higher order for this distribution cannot be derived easily, and may thus be approximated by Monte Carlo methods for practical applications.

Naturally, both extensions can also be combined by mixtures of shifted bivariate Skellam distributions. Parameter estimation for a shifted BSkellam distribution, a mixture of BSkellam distributions, and mixtures of shifted BSkellam distribution may be carried out by numerical maximization of the log likelihood. However, it may be noted that this procedure can become computationally demanding, because it requires optimization over a set of parameters with continuous and discrete-values entries.

Applications

In this application section, we present a Monte Carlo study in the next Section 4.1, followed by an application to real data in Section 4.2 and 4.3.

Monte Carlo study

In this section we briefly demonstrate some properties of the bivariate Skellam distribution and the parameter estimation procedures. More precisely, we consider two settings (termed Setting 1 and 2 in the following), dealing with a bivariate Skellam distribution with parameters (λ 0 , λ 1 , λ 2 ) = (4, 2, 3) and (λ 0 , λ 1 , λ 2 ) = (4, 35, 40), respectively. Figure 1 shows random samples with 500 observations each from Setting 1 in the left and Setting 2 in the right panel. The two settings are subject to different correlation, which equals 0.617 and 0.0966, respectively. For the two samples displayed, the estimated correlation equal 0.628 (p < 2.2e-16) and 0.0975 (p = 0.0292), respectively.

In a first Monte Carlo experiment, we investigate the departure from normality in the two settings. For this purpose, we generated 1000 samples having 1000 observations for each of the settings and performed the Shapiro-Wilk test for multivariate normality. In Setting 1, having comparably small values of λ i , i = 1, 2, 3, normality is rejected at 5%-level for 87.2% of the samples. However, the corresponding proportion in the second setting reduces to 30.1%. This indicates that the the bivariate Skellam can be approximated by a Gaussian distribution for large values of λ i , i = 1, 2, 3, which supports the results of Lemma 2.

In a second Monte Carlo experiment, we investigate the properties of the method of moments and maximum likelihood estimators. To this end, the two parameter sets described above were used for generating 1000 samples of observations of size n = 50, 100, and 200, respectively. Subsequently, the parameters were estimated by the method of moments and by maximization of the likelihood. Figure 2 shows Box plots of the estimated parameters. The method of moments estimates have white background and the boxes corresponding to the maximum likelihood estimator are coloured grey. The dashed black lines represent the true parameter values. Several phenomena are visible: First, both estimators seem to be consistent. Second, both estimators seem unbiased. Third, the maximum likelihood estimator has a lower variability than the method of moments estimator. Finally, it may be noted that normality of the estimators cannot be rejected for the large majority of samples at 5%-level (30 of 36, Shapiro-Wilk).

Soccer results

The data analyzed in this section are soccer results from the qualifications matches to the soccer world championship 2014 that took place in June and September 2012. The data are available for download on the FIFA web site (http://www.fifa.com/worldcup/preliminaries/matches/date=092012.html). More precisely, we consider the difference of goals scored by the home and the guest team during the first and second half, respectively, of the match. Figure 3 shows the 239 observations. The average difference in the first half equals 0.155 (s.d.: 1.19) and 0.172 (s.d: 1.54) in the second half. Moreover, the estimated correlation of the differences calculated from the first and second half takes the value 0.412 (p < 0.001). This indicates that an independent modelling of the variables would not be appropriate.

We fitted four different models to the data: first, a bivariate Skellam and a shifted bivariate Skellam distribution. Secondly, mixtures of bivariate Skellam and shifted bivariate Skellam distributions. These four models are denoted BSk, sBSk, mix-BSk, and mix-sBSk, respectively, in the following. The parameters of all models were estimated by maximization of the log-likelhood. Table 1 displays the maximized values of the log-likelihood and the model selection criteria AIC and BIC. Both mixture models have two components, because models with a higher number of components are rejected by both selection criteria.

The AIC shows a preference for the mix-sBSk, while the BIC prefers the more parsimonious mix-BSk without shift. The estimated parameters of the mix-BSk take the values ( λ01 , λ11 , λ21 ) = (0.383, 0.561, 0.955), ( λ02 , λ12 , λ22 ) = (2.190, 2.300, 1.490), (p 1 , p2 ) = (0.685, 0.315), and the respective estimates for mix-sBSk equal ( λ01 , λ11 , λ21 ) = (0.272, 0.463, 0.836), ( λ02 , λ12 , λ22 ) = (1.000, 1.120, 2.780), (p 1 , p2 ) = (0.503, 0.497), ( k11 , k21 ) = (0, -2), ( k12 , k22 ) = (0, 0). The simpler models with only one component are less favoured, which is underlined by the statistics presented in Table 2. This table shows the first three moments of the marginal distributions and the correlation for both the sample data and the estimated models.

The marginal means reproduced well by all models. As to standard deviation, skewness and correlation values, the mixture models show their superiority to the simpler models with only one component. Each of the two mixture models has its small advantages and disadvantages: While the model without shift comes closer to the correlation and skewness values of the data, the model incorporation a shift better captures the standard deviation of the sample.

The better fit of the mixture models also becomes visible in Figure 4, which shows histograms of the marginal distribution of the data with the four fitted models. From a contextual point of view, the application of a mixture model may be motivated by the comparably large heterogeneity of the matches carried out for a the qualification of the soccer world championship. For example, matches of more or less equally strong teams can be expected to have a distribution different to that of a matches of teams with large performance differences. However, further investigation of such hypothesis goes beyond the scope of this paper.

Analysis of EURO-Bund and Schatz future

In this section, we analyze data from the EURO-Bund and Schatz future, respectively, recorded on November 29 th 2011. The maturity of both futures was December 2011, and the smallest possible price change (tick) equals 0.01 and 0.005 for Bund and Schatz future, respectively. Prices of the two futures are recorded for intervals of two minute length, where the last trade in each interval determines the corresponding price. If no trades took place in an interval, the price was supposed to be unchanged. In order to avoid non-representative outliers, we restrict our observation period to the main trading hours from 9h15 to 17h15, resulting in 240 observations. During this period, the Bund future dropped slightly from 134.28 to 133.4, whereas the Schatz future remained almost unchanged (109.98 vs. 109.985). Figure 5 displays the price changes of the two futures in ticks, as well as their sample ACF, cross-correlation function, and a bivariate scatter plot. The two data series are time series, however, neither the ACF nor the cross-correlation function indicate the presence of any temporally lagged dependencies. Nevertheless, the two series are correlated at lag zero, the estimated correlation equals 0.664 (p < 0.001).

As in the previous sections, we fitted different models to the data. Here, we restrict ourselves to the simple bivariate Skellam (BSk) model and mixtures of bivariate Skellam distributions (mix-BSk), because preliminary analyses showed that a substantial increase in the log-likelihood could be easily achieved by increasing the number of mixture components, whereas shifted distributions only had a comparably small effect. Table 3 displays the maximized values of the log-likelihood and the model selection criteria AIC and BIC. Both AIC an BIC prefer the mixture model with 4 states.

The better fit of the mixture models also becomes visible in Figure 6, which shows histograms of the marginal distribution of the data with the two fitted models. The simple BSk fails to reproduce in particular the marginal of the Schatz future, for which the mix-BSk provides a more satisfactory image.

The estimated parameters of BSk take the values ( λ0 , λ1 , λ2 ) = (6.051, 5.672, 6.051), and the respective estimates for mix-BSk equal ( λ01 , λ11 , λ21 ) = (0.570, 1.375, 0.651), ( λ02 , λ12 , λ22 ) = (3.042, 0.000, 2.657), ( λ03 , λ13 , λ23 ) = (0.000, 5.873, 1.300), ( λ04 , λ14 , λ24 ) = (7.406, 0.000, 6.106),

(p 1 , p2 , p3 , p4 ) = (0.410, 0.283, 0.183, 0.125).

Proofs

Proof of Lemma 1.

-For any i ∈ {1, 2}, we have E(X i ) = λ i -λ 0 and V(X i ) = λ i + λ 0 . We have C ov (X 1 , X 2 ) = λ 0 . From these equalities follows the covariance matrix Σ. -Using the first point, we obtain cor(X 1 , X

2 ) = C ov (X 1 , X 2 )/ V(X 1 )V(X 2 ) = λ 0 / (λ 1 + λ 0 )(λ 2 + λ 0 ). -For any (t 1 , t 2 ) ∈ Z 2 , the characteristic function of (X 1 , X 2 ) is, φ (X1,X2) (t 1 , t 2 ) = E(e it1X1+it2X2 ) = E(e it1(U1-U0)+it2(U2-U0) ) = E(e it1U1+it2U2-i(t1+t2)U0 ) = φ U1 (t 1 )φ U2 (t 2 )φ U0 (-t 1 -t 2 )
= e λ1(e it 1 -1) e λ2(e it 2 -1) e λ0(e -it 1 -it 2 -1) .

This ends the proof of Lemma 1.

Proof of Lemma 2. Set W 1 = 1/(λ 1 -λ 0 )X 1 and W 2 = 1/(λ 2 -λ 0 )X 2 For λ 2 and λ 1 large enough, using the third point of Lemma 1 and

e x ∼ 0 1+x+x 2 /2 follows φ (W1,W2) (t 1 , t 2 ) = φ (X1,X2) (t 1 /(λ 1 -λ 0 ), t 2 /(λ 2 -λ 0 ))
= e λ1(e it 1 /(λ 1 -λ 0 ) -1) e λ2(e it 2 /(λ 2 -λ 0 ) -1) e λ0(e -it 1 /(λ 1 -λ 0 )-it 2 /(λ 2 -λ 0 ) -1)

≈ e λ1(it1/(λ1-λ0)-(1/2)t 2 1 /(λ1-λ0) 2 ) e λ2(it2/(λ2-λ0)-(1/2)t 2 2 /(λ2-λ0) 2 ) e λ0(-it1/(λ1-λ0)-it2/(λ2-λ0)-(1/2)(t1/(λ1-λ0)+t2/(λ2-λ0)) 2 ) = e it1 e it2 e -(1/2)(t 2 1 /(λ1-λ0))+t 2 2 /(λ2-λ0)+2t1t2/(λ1-λ0)(λ2-λ0))
for any (t 1 , t 2 ) ∈ Z 2 . The last term is the characteristic function of a bivariate random variable

Z ∼ N 2 1 1 , 1/(λ 1 -λ 0 ) λ 0 /(λ 1 -λ 0 )(λ 2 -λ 0 ) λ 0 /(λ 1 -λ 0 )(λ 2 -λ 0 ) 1/(λ 2 -λ 0 ) .
Therefore (W 1 , W 2 ) can be approximated by this distribution and, a fortiori,

(X 1 , X 2 ) ≈ N 2 λ 1 -λ 0 λ 2 -λ 0 , λ 1 -λ 0 λ 0 λ 0 λ 2 -λ 0 .
Lemma 2 is proved.

Proof of Lemma 3. For any (x, y) ∈ Z 2 holds

∂G(λ, z) ∂λ 1 = λ x 1 λ y 2 ∞ i=0 (λ 0 λ 1 λ 2 ) i (x + i)!(y + i)!i! = xλ x-1 1 λ y 2 ∞ i=0 (λ 0 λ 1 λ 2 ) i (x + i)!(y + i)!i! + λ x 1 λ y 2 ∞ i=1 i (λ 0 λ 2 ) i λ i-1 1 (x + i)!(y + i)!i! = x λ 1 G(λ, z) + λ x 1 λ y 2 ∞ i=0 (λ 0 λ 2 ) i+1 λ i 1 (x + i + 1)!(y + i + 1)!i! = x λ 1 G(λ, z) + λ 0 λ 1 G(λ, (x + 1, y + 1)).
The other equalities can be proved similarly, and the proof of Lemma 3 is complete.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -8 -4 0 2 4 6 -5 0 5

x 1

x 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 10 20 30 40 50 20 30 40 50

x 1

x 2 Fig. 1 Two samples with 500 observations each, generated from a bivariate Skellam distribution with true parameters (λ 0 , λ 1 , λ 2 ) = (4, 2, 3) and (λ 0 , λ 1 , λ 2 ) = (4, 35, 40) in the left (respectively right) panel. To improve the visualization, Gaussian jitter was added to the observations. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 2 3 4 5 6 7 8 λ 0

Setting 1 n = 50 100 200 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 2 4 6 λ 1 50 100 200 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 2 4 6 λ 2 50 100 200 M.M. est. M.L. est. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -20 -10 0 10 20 λ 0

Setting 2 n = 50 100 200 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 10 20 30 40 50 λ 1 50 100 200 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 20 30 40 50 60 λ 2 50 100 200

Fig. 2 Parameter estimates from 100 from simulated sequences of length n = 50, 100, and 200. The upper three panels were simulated from a bivariate Skellam with parameters (λ 0 , λ 1 , λ 2 ) = (4, 2, 3), the lower three panels result from the parameters (λ 0 , λ 1 , λ 2 ) = (4, 35, 40) in the left and right panel, respectively. Box plots corresponding to the method of moments (M.M.) and maximum likelihood (M.L.) estimator have white (respectively grey) background. The dashed black lines represent the true parameter values.

Fig. 3 Goal difference between home and guest team at half-time and at the end of the match. For improved visualization, Gaussian jitter was added to the observations. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -4 -2 0 2 4 6 -4 -2 0 2 4 Goal difference 1st half Goal difference 2nd half q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 2nd half q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q BSk sBSk mix-BSk mix-sBSk Series dat$t2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -15 -10 -5 0 5 10

Fig. 4

 4 Fig. 4 Histogram of the marginal distributions with fitted models. 1st half

Fig. 5

 5 Fig.5Price changes of Bund and Schatz future in the top panels. The middle panels show the sample ACF of the two data series. The lower two panels show a bivariate scatter plot of the data and the empirical cross-correlation. For improved visualization, Gaussian jitter was added to the observations in the scatter plot.

Fig. 6

 6 Fig.6Histogram of the marginal distributions with fitted models.Bund future

Table 1

 1 Model selectionThis table shows log-likelihood, AIC, and BIC values for the four estimated models.

		BSk sBSk mix-BSk mix-sBSk
	logL	-814	-802	-794	-787
	AIC	1634	1614	1600	1595
	BIC	1645	1632	1621	1630

Table 2

 2 Moments and correlationThis table summarizes the first three moments of the marginal distributions and the correlation of the data and models.

		Mean	S.d.	Skew.		Corr.
	Data	0.155 0.172 1.19 1.54	0.00101	-0.439	0.412
	BSk	0.155 0.172	2.1 2.12	0.0507 0.0555	0.461
	sBSk	0.155 0.172 1.49 2.51	0.0845	0.292	0.346
	mix-BSk	0.155 0.172 1.43 1.56 -0.00228	-0.503	0.434
	mix-sBSk 0.155 0.172 1.19 1.61	0.0422	-0.102	0.338

Table 3

 3 Model selectionThis table shows log-likelihood, AIC, and BIC values for the four estimated models.

	no. states	1	2	3	4	5
	logL	-1075 -1149 -1096 -1075 -1073
	AIC	2173	2310	2211	2173	2176
	BIC	2215	2331	2242	2215	2228

Acknowledgements

We would like to thank the Chair of Quantitative Finance of the Ecole Centrale Paris for granting us access to their data base. Special thanks go to F. Abergel for his suggestion to closer investigating futures on government bonds and R. Zaatour for his continuous support with the data preparation process. Furthermore, we are would like to render thanks to two anonymous reviewers for their helpful comments.