Joachim Lebovits
email: jolebovits@gmail.com.

Mark Podolskij
email: mpodolskij@math.au.dk

Estimation of the global regularity of a multifractional Brownian motion

Keywords: consistency, Hurst parameter, multifractional Brownian motion, power variation AMS Subject Classification: 60G15, 60G22 62G05, 62M09, 60G17

This paper presents a new estimator of the global regularity index of a multifractional Brownian motion. Our estimation method is based upon a ratio statistic, which compares the realized global quadratic variation of a multifractional Brownian motion at two different frequencies. We show that a logarithmic transformation of this statistic converges in probability to the minimum of the Hurst function, which is, under weak assumptions, identical to the global regularity index of the path.

Introduction

Fractional Brownian motion (fBm) is one of the most prominent Gaussian processes in the probabilistic and statistical literature. Popularized by Mandelbrot and van Ness [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF] in 1968, it found various applications in modeling stochastic phenomena in physics, biology, telecommunication and finance among many other fields. Fractional Brownian motion is characterized by its self-similarity property, the stationarity of its increments and by its ability to match any prescribed constant local regularity. Mathematically speaking, for any H ∈ (0, 1), a fBm with Hurst index H, denoted by B H = (B H t ) t≥0 , is a zero mean Gaussian process with the covariance function given by

E[B H s B H t ] = 1 2 t 2H + s 2H -|t -s| 2H .
Various representations of fBm can be found in the existing literature; we refer to [START_REF] Nualart | The Malliavin Calculus and Related Topics[END_REF][START_REF] Lebovits | Stochastic integration with respect to multifractional Brownian motion via tangent fractional Brownian motions[END_REF] and references therein. The Hurst parameter H ∈ (0, 1) determines the path properties of the fBm: (i) The process (B H t ) t≥0 is self-similar with index H, i.e. (a H B H t ) t≥0 = (B H at ) t≥0 in distribution, (ii) (B H t ) t≥0 has Hölder continuous paths of any order strictly smaller than H, (iii) fractional Brownian motion has short memory if an only if H ∈ (0, 1/2]. Moreover, fBm presents long range dependance if H belongs to (1/2, 1). The statistical estimation of the Hurst parameter H in the high frequency setting, i.e. the setting of mesh converging to 0 while the Definition 1. Let h : R + → (0, 1) be a deterministic function and B be a fractional Brownian field. A multifractional Brownian motion (mBm) with functional parameter h is the Gaussian process B h = (B h t ) t∈R + defined by B h t := B(t, h(t)), for all t ∈ R + .

Define, for any x in (0, 1), the positive real c x by setting:

c x := 2 cos(πx)Γ(2 -2x) x(1 -2x) 1 2 , (1.1)
where Γ denotes the standard gamma function. For any function h : R + → (0, 1), it is easy to verify that the process B h := (B h t ) t∈R + defined by

B h t = 1 c h(t) R exp(itx) -1 |x| h(t)+1/2 W (dx), (1.2)
where W denotes a complex Gaussian measure1 , is a multifractional Brownian motion with functional parameter h. Intuitively speaking, the multifractional Brownian motion behaves locally as fractional Brownian motion, but the functional parameter h is time-varying. Moreover, it remains linked to local regularity of B h , but in a less simple way than in the case of the fBm. More precisely, if we assume that h belongs to the set C η ([0, 1], R), for some η > 0, and is such that

0 < h min := min t∈[0,1] h(t) ≤ h max := max t∈[0,1] h(t) < min{1, η}, (1.3)
then h min is the regularity parameter of B h (see [ACLV00, Corollaries 1,2 and Proposition 10]).

In this setting the functional parameter h needs to be estimated locally in order to get a full understanding of the path properties of the multifractional Brownian motion B h . Bardet and Surgailis [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional Gaussian processes[END_REF] have proposed to use a local power variation of higher order filters of increments of B h to estimate the function h. More specifically, they prove the law of large numbers and a central limit theorem for the local estimator of h (i) based on log-regression of the local quadratic variation, (ii) based on a ratio of local quadratic variations.

In this paper we are aiming at the estimation of the parameter h min , which represents the regularity (or smoothness) of the multifractional Brownian motion B h = (B h t ) t≥0 . For this particular statistical problem the local estimation approach investigated in [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional Gaussian processes[END_REF] appears to be rather inconvenient. Instead our method relies on a ratio statistic, which compares the global quadratic variation at two different frequencies. We remark that in general it is impossible to find a global rate a n such that the normalized power variation a n V (B h , p) n T converges to a nontrivial limit. However, ratios of global power variations can very well be useful for statistical inference. Indeed, we will show that under appropriate conditions on the functional parameter h, the convergence

S n (B h ) := n-1 i=0 B h i+1 n -B h i n 2 n-2 i=0 B h i+2 n -B h i n 2 -→ n→+∞ 2 -2h min ,
holds in probability.

Then a simple log transformation gives a consistent estimator of the global regularity h min of a mBm.

The paper is structured as follows. Section 2 presents the basic distribution properties of the multifractional Brownian motion, reviews the estimation methods from [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional Gaussian processes[END_REF] and states the main asymptotic results of the paper. Proofs are given in Section 3.

Background and main results

In [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional Gaussian processes[END_REF] Bardet and Surgailis deal with a little bit more general processes than multifractional Brownian motions. However, in order not to overload the notations we will focus in this paper on the normalized multifractional Brownian motion (i.e. the mBm defined by (1.2)). From now on we will refer to this process as the multifractional Brownian motion and denote it by

B h = (B h t ) t≥0 .

Basic properties and local estimation of the functional parameter h

We start with the basic properties of the mBm B h with functional parameter h. Its covariance function is given by the expression

R h (t, s) := E[B h t B h s ] = c 2 ht,s 2c h(t) c h(s) |t| 2ht,s + |s| 2ht,s -|t -s| 2ht,s , (2.1)
where h t,s := h(t)+h(s) 2 and c x has been defined in (1.1). It is easy to check that x → c x is a C ∞ ((0, 1))-function. The local behaviour of the multifractional Brownian motion is best understood via the relationship

u -h(t) (B h t+us -B h t ) s≥0 f.d.d. -→ B h(t) s s≥0
as u → 0, where f.d.d.

-→ denotes the convergence of finite dimensional distributions. Hence, in the neighbourhood of any t in (0, 1), the mBm B h behaves as fBm with Hurst parameter h(t). This observation is essential for the local estimation of the functional parameter h. In the following we will briefly review the statistical methods of local inference investigated in Bardet and Surgailis [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional Gaussian processes[END_REF], which is based on high frequency observations B h 0 , B h 1/n , . . . , B h (n-1)/n , B h 1 . While the original paper is investigating rather general Gaussian models whose tangent process is a fractional Brownian motion, we will specialize their asymptotic results to the framework of multifractional Brownian motion.

Let us introduce the generalized increments of a process Y = (Y t ) t≥0 . Consider a vector of coefficients a = (a 0 , . . . , a q ) ∈ R q+1 and a natural number m ≥ 1 such that q j=0 j k a j = 0 for k = 0, . . . , m -1 and q j=0 j m a j = 0.

In this case the vector a ∈ R q+1 is called a filter of order m. The generalised increments of Y associated with filter a at stage i/n are defined as

∆ n i,a Y := q j=0 a j Y i+j n .
Standard examples are a (1) = (-1, 1), ∆ n i,a (1) Y = Y (i+1)/n -Y i/n (first order differences) and

a (2) = (1, -2, 1), ∆ n i,a (2) Y = Y (i+2)/n -2Y (i+1)/n + 2Y i/n
(second order differences). In both cases we have that q = m. Now, we set ψ(x, y) := (|x + y|)/(|x| + |y|) and set

Λ(H) := E[ψ(∆ n 0,a B H , ∆ n 1,a B H )],
H ∈ (0, 1).

The function Λ does not depend on n and is strictly increasing on the interval (0, 1). For any α ∈ (0, 1), which determines the local bandwidth, the ratio type estimator of h(t) is defined as

h n,α t := Λ -1 1 card{i ∈ 0, n -q -1 : |i/n -t| ≤ n -α } i∈ 0,n-q-1 : |i/n-t|≤n -α ψ(∆ n i,a B h , ∆ n i+1,a B h ) .
(2.2) Here and throughout the paper we denote p, q := {p, p + 1, p + 2, . . . , q} for any p, q ∈ N with p ≤ q. The authors of [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional Gaussian processes[END_REF] only investigate the estimator h n,α t relative to the filter a = a (2) , which we assume in this subsection from now on. The consistency and asymptotic normality of the estimator h n,α t is summarized in the following theorem. We remark that the condition for the central limit theorem crucially depends on the interplay between the bandwidth parameter α and the Hölder index η of the function h.

Theorem 2.1. ([BS13, Proposition 3]) Assume that h belongs to

C η ([0, 1]) and that condition (1.3) is satisfied. (i) For any t ∈ (0, 1) and α ∈ (0, 1) it holds that h n,α t P -→ h(t), as n → ∞. (ii) When α > max 1 1+2 min(η,2) , 1 -4(min(η, 2) -sup t∈(0,1) h(t)) it holds that √ 2n 1-α h n,α t -h(t) d -→ N (0, τ 2 ) as n → ∞,
where the asymptotic variance τ 2 is defined in [BS13, Eq. (2.17)].

The paper [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional Gaussian processes[END_REF] contains the asymptotic theory for a variety of other local estimators of h(t). We dispense with the detailed exposition of these estimators, since only h n,α t is somewhat related to our estimation method.

Remark 1. Nowadays, it is a standard procedure to consider higher order filters for Gaussian processes to obtain a central limit theorem for the whole range of Hurst parameters. Let us shortly recall some classical asymptotic results, which are usually referred to as Breuer-Major central limit theorems. We consider the scaled power variation of a fractional Brownian motion B H with Hurst parameter H ∈ (0, 1) based on first order filter a (1) and second order filter a (2) :

V (B H , p; a (1) ) n := n -1+pH n-1 i=0 |∆ n i,a (1) B H | p and V (B H , p; a (2) ) n := n -1+pH n-2 i=0 |∆ n i,a (2) B H | p .
It is well known that, after an appropriate normalization, the statistic V (B H , p; a (1) ) n exhibits asymptotic normality for H ∈ (0, 3/4], while it converges to the Rosenblatt distribution for H ∈ (3/4, 1). On the other hand, the statistic V (B H , p; a (2) ) n exhibits asymptotic normality for all H ∈ (0, 1). We refer to [START_REF] Breuer | Central limit theorems for nonlinear functionals of gaussian fields[END_REF][START_REF] Taqqu | Convergence of integrated processes of arbitrary hermite rank[END_REF] for a detailed exposition.

Estimation of the global regularity parameter h min

In this section we will construct a consistent estimator of the global regularity parameter h min , which has been defined at (1.3). Our first condition is on the set h -1 ({h min }), which is necessarily compact since h belongs to C η ([0, 1]). We assume that this set has the following form

M h := h -1 ({h min }) = q i=1 [a i , b i ]   m j=1 {x j }   , (q, m) ∈ N 2 \ (0, 0), (2.3)
where N = {0, 1, 2, . . .} and the intervals [a i , b i ] are disjoint and such that none of the x j 's belongs to

q i=1 [a i , b i ].
Depending on whether q ≥ 1 or q = 0, we will need an additional assumption. Below, we denote by h (A ) There exist positive integers p j such that function h is p j times continuously left and right differentiable at point x j for j = 1, . . . , m such that

p j = min{p : h (p) l (x j ) = 0} = min{p : h (p) r (x j ) = 0}.
We remark that since h reaches its minimum at points x j , we necessarily have that h

(p j ) r (x j ) > 0 and that h (p j ) l (x j ) > 0 if p j is even and h (p j ) l (x j ) < 0 if p is odd. Now,

we proceed with the construction of the consistent estimator of the global regularity parameter h min based on high frequency observations

B h 0 , B h 1/n , . . . , B h (n-1)/n , B h 1 .
First of all, let us remark that considering the estimator min t∈[0,1] h n,α t , where h n,α t has been introduced in the previous section, is not a trivial matter since the functional version of Theorem 2.1 is not available. Instead our statistics relies on the global quadratic variation rather than local estimates.

For the mBm B h = (B h t ) t∈[0,1] , we introduce the notations

V (B h ; k) n := n-k i=0 B h i+k n -B h i n 2 , S n (B h ) := V (B h ; 1) n V (B h ; 2) n .
(2.4) Our first result determines the limit of

E[V (B h ; 1) n ]/E[V (B h ; 2) n ]. Proposition 2.2. Let h : [0, 1] → (0, 1) be a deterministic C η ([0, 1])-function satisfying (1.
3) and such that the set M h has the form (2.3). If q = 0 we also assume that condition (A ) holds. Define

U h n := E[V (B h ; 1) n ] E[V (B h ; 2) n ] .
Then it holds that

lim n→+∞ U h n = 1 2 2h min . (2.5)
The convergence result of Proposition 2.2 is rather intuitive when q ≥ 1, which means that the minimum of the function h is reached on a set of positive Lebesgue measure. In this setting it is quite obvious that the statistic

V (B h ; k) n is dominated by squared increments (B h (i+k)/n -B h i/n ) 2 for i/n ∈ ∪ q i=1 [a i , b i ].
Thus, the estimation problem is similar to the estimation of the Hurst parameter of a fractional Brownian motion (B h min t

) t∈∪ q i=1 [a i ,b i ]
with Hurst parameter h min , for which the convergence at (2.5) is well known. When q = 0, and hence Leb(M h ) = 0, the proof of Proposition 2.2 becomes much more delicate.

Remark 2. Assume for illustration purpose that q = 0, m = 1, x := x 1 and p := p 1 . Condition (A ) is crucial to determine the precise asymptotic expansion of the quantity E[V (B h ; k) n ]. The lower and upper bounds in (3.16) and (3.17) in the proof show that

E[V (B h ; k) n ] = O n 1-2h min (ln n) 1/p as n → +∞, for k = 1, 2.
The condition min{p :

h (p) l (x) = 0} = min{p : h (p)
r (x) = 0} of assumption (A ) is not essential for the proofs. For instance, when min{p :

h (p) l (x) = 0} > min{p : h (p) r (x) = 0} the expectation E[V (B h ; k) n ]
would be dominated by the terms in the small neighbourhood on the right hand side of x and the statement of Proposition 2.2 can be proved in the same manner.

Our main result shows that the statistic S n (B h ) and the quantity U h n are asymptotically equivalent in probability.

Theorem 2.3. Assume that h ∈ C 2 ([0, 1]) and the set M h has the form (2.3). If q = 0 we also assume that condition (A ) holds. Then we have the following result:

S n (B h ) P -→ 1 2 2h min . (2.6)
In particular, the following convergence holds:

h min := - ln(S n (B h )) 2 ln(2) P -→ h min .
(2.7)

The asymptotic result of Theorem 2.3 can be extended to more general Gaussian processes than the mere multifractional Brownian motion. As it has been discussed in [START_REF] Bardet | Nonparametric estimation of the local Hurst function of multifractional Gaussian processes[END_REF], when a Gaussian process possesses a tangent process B h(t) at time t, we may expect Theorem 2.3 to hold under certain assumptions on its covariance kernel. We refer to assumptions (A) κ and (B) α therein for more details on sufficient conditions.

The rate of convergence, or a weak limit theorem, associated with the consistency result at (2.7) is a more delicate issue. In the setting q = 0, which implies that Leb(M h ) = 0, the bias associated with the convergence in (2.5) may very well dominate the variance of the estimator. A careful inspection of our proof, and more specifically of statement (3.8) and Remark 3, implies that the bias in Proposition 2.2 has a logarithmic rate. Thus, weak limit theorems for the estimator h min are out of reach in this framework. When q ≥ 1 and hence Leb(M h ) > 0, one may hope to find better rates of convergence for the estimator h min . However, we dispense with the exact exposition of this statistical problem.

Proofs

Throughout this section we denote all positive constants by C, or C p if they depend on an external parameter p, although they may change from line to line.

Proof of Proposition 2.2

For k = 1, 2 we introduce the notation

V (k) n := n-k i=0 k n 2h(i/n) , (3.1)
which serves as the first order approximation of the quantity

E[V (B h ; k) n ]. Applying [BS10, Lemma 1 p.13] we conclude that E[V (B h ; k) n ] -V (k) n ≤ C ln n n η∧1 n-k i=0 i n 2h(k/n) ≤ C ln n n 2h min -1+η∧1 (3.2)
for any (n, k) ∈ N × {1, 2}. We have the inequality

U h n - 1 2 2h min ≤ |E[V (B h ; 1) n ] -V (1) n | + |E[V (B h ; 2) n ] -V (2) n | V (2) n + V (1) n V (2) n - 1 2 2h min =: µ (1) n + µ (2) n . (3.3)
We first show that µ

(1) n → 0 as n → ∞. When h min = h max we trivially have µ

(1) n = 0. If h min < h max , we fix ǫ ∈ (0, h maxh min ). By Leb(A) we denote the Lebesgue measure of any measurable set A. We have that Leb h -1 ([h min , h min + ǫ]) > 0.

Thus, there exists n 0 ∈ N such that for all n ≥ n 0 it holds that

Card{i ∈ 0, n -k ; h(i/n) ∈ [h min , h min + ǫ]} ≥ n Leb h -1 ([h min , h min + ǫ]) /2.
This implies that

V (2) n ≥ i∈ 0,n-k ; h(i/n)∈[h min ,h min +ǫ] 2 n 2h(i/n) ≥ Cn 1-2(h min +ǫ) .
Hence, applying Inequality (3.2), we conclude that:

µ (1) n ≤ C ln n • n 2ǫ-η∧1 ,
which proves that µ

(1) n → n→+∞ 0, for any ǫ small enough.

Convergence of µ

(2) n in the case q ≥ 1 We first prove that µ

(2) n → 0 in the case q ≥ 1. Assume again that h min < h max . First, we observe the lower bound

V (k) n ≥ q l=1 i∈ 0,n-k ; i/n∈[a l ,b l ] k n 2h(i/n) = k n 2h min q l=1 card{i ∈ 0, n -k ; i/n ∈ [a l , b l ]} ≥ n k n 2h min q l=1 b l -a l -2 n . (3.4)
For the upper bound we fix 0 < ǫ < h maxh min and consider the decomposition

V (k) n = i∈ 0,n-k ; h(i/n)∈[h min ,h min +ǫ] k n 2h(i/n) + i∈ 0,n-k ; h(i/n) ∈[h min ,h min +ǫ] k n 2h(i/n) . Setting λ n (ǫ) := n -1 card{i ∈ 0, n -k ; h(i/n) ∈ [h min , h min + ǫ]}, we deduce the assertions λ n (ǫ) → Leb h -1 ([h min , h min + ǫ]) as n → ∞, Leb h -1 ([h min , h min + ǫ]) → Leb h -1 ({h min }) = q l=1 (b l -a l ) > 0 as ǫ → 0.
Now, we conclude that

V (k) n ≤ nλ n (ǫ) k n 2h min + n(1 -λ n (ǫ)) k n 2(h min +ǫ) . (3.5)
Throughout the proofs we write lim for lim inf and lim for lim sup. Applying inequalities (3.4) and (3.5), we obtain that

lim n→+∞ n q l=1 (b l -a l -2 n ) nλ n (ǫ) + n(1 -λ n (ǫ)) 2 n 2ǫ ≤ lim n→+∞ 2 2h min V (1) n V (2) n ≤ lim n→+∞ 2 2h min V (1) n V (2) n ≤ lim n→+∞ nλ n (ǫ) + n(1 -λ n (ǫ)) 1 n 2ǫ n q l=1 (b l -a l -2 n )
.

Hence, we deduce that

2 -2h min Leb h -1 ({h min }) Leb (h -1 ([h min , h min + ǫ])) ≤ lim n→+∞ V (1) n V (2) n ≤ lim n→+∞ V (1) n V (2) n ≤ 2 -2h min Leb h -1 ([h min , h min + ǫ]) Leb (h -1 ({h min }))
.

By letting ǫ tend to 0, we readily deduce taht µ

(2) n → 0 as n → +∞.

Convergence of µ

(2) n in the case q = 0 Without loss of generality we assume that m = 1 and M h = h -1 ({h min }) = {x} with x ∈ (0, 1). Recall that in this setting we assume condition (A ) with p := p 1 . We let γ be a positive number such that γ < 2 -1 min{|h

(p) l (x)|, h (p)
r (x)}. Now, there exists a ǫ = ǫ(γ) > 0 with ǫ < min{x, 1x, γ} such that: ∀y > x with 0 < yx < ǫ :

h min + 1 p! (y -x) p (h (p) r (x) -γ) ≤ h(y) ≤ h min + 1 p! (y -x) p (h (p) r (x) + γ), (3.6)
∀y < x with 0 < xy < ǫ :

h min + 1 p! (y -x) p (h (p) l (x) -(-1) p γ) ≤ h(y) ≤ h min + 1 p! (y -x) p (h (p) l (x) + (-1) p γ). (3.7)
We proceed with the derivation of upper and lower bounds for the quantity µ

(2)

n . We start with the decomposition yε) , where we have set

V (k) n = Γ (1) n,k (γ, ǫ) + Γ (2) n,k (γ, ǫ) + Γ (3) n,k (γ, ǫ) where Γ (1) n,k (γ, ǫ) := i∈ 0,n-k ; i/n∈[x,x+ǫ] k n 2h(i/n) ; Γ (2) n,k (γ, ǫ) := i∈ 0,n-k ; i/n∈[x-ǫ,x) k n 2h(i/n) ; Γ (3) n,k (γ, ǫ) := i∈ 0,n-k ; i/n∈[x-ǫ,x+ǫ] c k n 2h(i/n) . It is clear that Γ (3) n,k (γ, ǫ) ≤ n(k/n) 2h(
y ǫ := argmin{h(u) : u ∈ (x -ǫ, x + ǫ) c ∩ [0, 1]}.
For the other two quantities, we deduce that Γ (r)

n,k (γ, ǫ) ≤ Γ (r) n,k (γ, ǫ) ≤ Γ (r) n,k (γ, ǫ) with Γ (1) n,k (γ, ǫ) := k n 2h min i∈ 0,n-k : i/n∈[x,x+ǫ] k n 2(p!) -1 (i/n-x) p (h (p) r (x)+γ) , Γ (2) n,k (γ, ǫ) := k n 2h min i∈ 0,n-k : i/n∈[x-ǫ,x) k n 2(p!) -1 (i/n-x) p (h (p) l (x)+(-1) p γ) and Γ (1) n,k (γ, ǫ) := Γ (1) n,k (-γ, ǫ) and Γ (2) n,k (γ, ǫ) := Γ (2)
n,k (-γ, ǫ). Using (3.6) and (3.7), it is easy to see that, for every (k, n) ∈ {1, 2} × N:

µ (2) n (γ, ǫ) ≤ V (1) n V (2) n ≤ µ (2) n (γ, ǫ), (3.8) with µ (2) n (γ, ǫ) := Γ (1) n,1 (γ, ǫ) + Γ (2) n,1 (γ, ǫ) Γ (1) n,2 (γ, ǫ) + Γ (2) n,2 (γ, ǫ) + n(2/n) 2h(yε) , µ (2) n (γ, ǫ) := Γ (1) n,1 (γ, ǫ) + Γ (2) n,1 (γ, ǫ) + n(1/n) 2h(yε) Γ (1) n,2 (γ, ǫ) + Γ (2) n,2 (γ, ǫ) . From (3.8) we obtain that 0 ≤ 2 2h min µ (2) n ≤ 2 2h min V (1) n V (2) n -1 ≤ U n (γ, ǫ) + U n (-γ, ǫ), ( 3.9)
where

U n (γ, ǫ) := |∆ n,2 (γ, ǫ)| -1 |2 2h min ∆ n,1 (γ, ǫ) -∆ n,2 (γ, ǫ)| + 2n 1-2h(yǫ) , (3.10) ∆ n,k (γ, ǫ) := Γ (1) n,k (γ, ǫ) + Γ (2) n,k (γ, ǫ), ∆ n,k (γ, ǫ) := ∆ n,k (-γ, ǫ). (3.11)
In view of (3.9) it is sufficient to show that lim

γ→0 lim n→+∞ U n (γ, ǫ) = 0. Define d γ := 2(p!) -1 (h (p) r (x) + γ) and d ′ γ := 2(p!) -1 (h (p) l (x) + (-1) p γ).
For any (a, b) in R + × (R \ {0}), we also set

S n,k (a, ǫ) := i∈ 0,n-k : i/n∈[x,x+ǫ] k n a(i/n-x) p , (3.12) T n,k (b, ǫ) := i∈ 0,n-k : i/n∈[x-ǫ,x) k n b(i/n-x) p . (3.13)
We deduce the identities Γ

(1)

n,k (γ, ǫ) = (k/n) 2h min S n,k (d γ , ǫ) and Γ (2) n,k (γ, ǫ) = (k/n) 2h min T n,k (d ′ γ , ǫ).

Note moreover that d ′

γ > 0 when p is even and d ′ γ < 0 when p is odd. We therefore assume from now on that b > 0 when p is even and that b < 0 when p is odd. For any η ∈ R \ {0}, we define

f (η) n,k (u) := k n η(u-x) p . Since i → f (a) n,k (i/n) is decreasing on [nx] + 1, [n(x + ǫ)] while i → f (b) n,k (i/n) is increasing if p even (resp. decreasing if p odd) on [n(x -ǫ)] + 1, [nx]
, one can use an integral test for convergence, which provides us with the following upper bounds n βn(a) αn(a) y 1/p-1 e -y dy p(a ln(n/k))

1/p ≤ S n,k (a, ǫ) ≤ n µn(a)
τn(a) y 1/p-1 e -y dy p(a ln(n/k)) 1/p , (3.14)

n β ′ n (b) α ′ n (b) y 1/p-1 e -y dy -ρ (b) n,k (ǫ) p((-1) p b ln(n/k)) 1/p ≤ T n,k (b, ǫ) ≤ n µ ′ n (b) τ ′ n (b) y 1/p-1 e -y dy -ρ (b) n,k (ǫ) p((-1) p b ln(n/k)) 1/p . (3.15)
Here we use the notation

α n (a) := a ln(n/k) [nx] + 1 n -x p , β n (a) := a ln(n/k) [n(x + ǫ)] + 1 n -x p , τ n (a) := a ln(n/k) [nx] n -x p , µ n (a) := a ln(n/k) [n(x + ǫ)] n -x p and ρ (b) n,k (ǫ) := f (b) n,k ( [n(x-ǫ)]+1 n ) + f (b) n,k ( [nx] n ). Furthermore, (α ′ n (b), β ′ n (b), τ ′ n (b), µ ′ n (b)) := (z (1) n (b), z (2) n (b), z (3) n (b), z (4) n (b)) if p is even, (α ′ n (b), β ′ n (b), τ ′ n (b), µ ′ n (b)) := (z (3) n (b), z (4) n (b), z (1) n (b), z (2) n (b)) if p is odd,
where we have set

z (1) n (b) := b ln(n/k) [nx] -2 n -x p , z (2) n (b) := b ln(n/k) [n(x -ǫ)] + 1 n -x p , z (3) n (b) := b ln(n/k) [nx] -1 n -x p , z (4) n (b) := b ln(n/k) [n(x -ǫ)] + 2 n -x p .
In view of the inequalities (3.14) and (3.15), as well as identities (3.12) and (3.13), we then deduce that

n 1-2h min k 2h min u n,k,p (d γ ) (ln(n/k)) 1/p • 1 d γ ≤ Γ (1) n,k (γ, ǫ) ≤ n 1-2h min k 2h min v n,k,p (d γ ) (ln(n/k)) 1/p • 1 d γ , (3.16) n 1-2h min k 2h min u ′ n,k,p (d ′ γ ) (ln(n/k)) 1/p • 1 |d ′ γ | ≤ Γ (2) n,k (γ, ǫ) ≤ n 1-2h min k 2h min v ′ n,k,p (d ′ γ ) (ln(n/k)) 1/p • 1 |d ′ γ | . (3.17)
Here we have used the notation

y 1/p-1 e -y dy - (-1) p b ln(n/k) 1/p ρ (b) n,k (ǫ) pn , v ′ n,k,p (b) := 1 p µ ′ n (b) τ ′ n (b) y 1/p-1 e -y dy - (-1) p b ln(n/k) 1/p ρ (b) n,k (ǫ) pn . Since Γ (r) n,k (γ, ǫ) = Γ (r)
n,k (-γ, ǫ), (3.16) and (3.17) also provide us with upper and lower bounds for Γ (r) n,k (γ, ǫ). Finally, we obtain the following lower and upper bounds

n 1-2h min k 2h min (ln(n/k)) 1/p • Λ n,k (γ, ǫ) ≤ ∆ n,k (γ, ǫ) ≤ n 1-2h min k 2h min (ln(n/k)) 1/p Λ ′ n,k (γ, ǫ), (3.18) n 1-2h min k 2h min (ln(n/k)) 1/p • Λ n,k (-γ, ǫ) ≤ ∆ n,k (γ, ǫ) ≤ n 1-2h min k 2h min (ln(n/k)) 1/p Λ ′ n,k (-γ, ǫ), (3.19) where Λ n,k (γ, ǫ) := 1 d γ • u n,k,p (d γ ) + 1 |d ′ γ | • u ′ n,k,p (d ′ γ ), Λ ′ n,k (γ, ǫ) := 1 d γ • v n,k,p (d γ ) + 1 |d ′ γ | • v ′ n,k,p (d ′ γ ).
Denote c p := +∞ 0 y 1/p-1 e -y dy. Recalling the definition of the constants d γ and d ′ γ , a straightforward computation shows that, for any (k,

k ′ ) ∈ {1, 2} 2 with k = k ′ : lim n→+∞ Λ n,k (γ, ǫ) = lim n→+∞ Λ ′ n,k (γ, ǫ) = c p p (1/d γ + 1/|d ′ γ |), (3.20) lim n→+∞ |Λ ′ n,k ′ (γ, ǫ) -Λ n,k (-γ, ǫ)| ≤ C (2|γ| + |1/d -γ -1/d γ + 1/|d ′ -γ | -1/|d ′ γ ||) ≤ C|γ|. (3.21)
Starting from (3.18), and using (3.20) and (3.21), we see that there exists a positive integer n 0 and C > 0 such that for all n ≥ n 0

|∆ n,k (γ, ǫ)| -1 ≤ C (ln(n/k)) 1/p n 1-2h min k 2h min . (3.22)
Finally, inequalities (3.20), (3.21) and (3.22) imply that there exists a positive integer N such that for all n ≥ N :

U n (γ, ǫ) ≤ C |γ| + (ln n) 1/p n 2(h(yǫ)-h min ) .
From the previous inequality, lim n→+∞ U n (γ, ǫ) ≤ C|γ| and thus we get lim γ→0 lim n→+∞ U n (γ, ǫ) = 0, which completes the proof. Remark 3. In the previous proof (in the case q = 0), using (3.20), one can also see that the bias related to the convergence of µ

(2) n to 0 is of order 1/ ln n.

Proof of Theorem 2.3

In the first step we will find an upper bound for the covariance function of the increments of B h . We define

r n (i, j) := cov B h i+k n -B h i n , B h j+k n -B h j n , k = 1, 2.
Recalling the notation at (2.1), we conclude the identity

r n (i, j) = R h i + k n , j + k n -R h i n , j + k n -R h i + k n , j n + R h i n , j n .
Since h ∈ C 2 ([0, 1]) and the function c defined at (1.1) is a C ∞ ((0, 1))-function, we deduce by an application of Taylor expansion

|r n (i, j)| ≤ n -2 2 l,l ′ =1 |∂ ll ′ R h (ψ n ij )| for |i -j| > 2, (3.23)
where ∂ ll ′ R h denotes the second order derivative in the direction of x l and x l ′ , and ψ n ij ∈ (i/n, (i + k)/n) × (j/n, (j + k)/n). Now, we will compute an upper bound for the right side of (3.23) for i = j. First, we observe that

R h (t, s) = F (t, s) G(t, s, h(t) + h(s)),
where

F (t, s) = c 2 ht,s c h(t) c h(s) , G(t, s, H) = 1 2 |t| H + |s| H -|t -s| H .
We remark that G(t, s, 2H) is the covariance kernel of the fractional Brownian motion with Hurst parameter H ∈ (0, 1). Since h ∈ C 2 ([0, 1]), c ∈ C ∞ ((0, 1)) and c x = 0 for x ∈ (0, 1), we conclude that

|∂ l F (t, s)|, |∂ ll ′ F (t, s)| ≤ C, l, l ′ = 1, 2, (t, s) ∈ [0, 1] 2 .
We concentrate on the second order derivative ∂ 11 R h (ψ n ij ); the estimates for the other second order derivatives are obtained similarly. We have that

∂ 11 R h (t, s) = ∂ 11 F (t, s) • G(t, s, h(t) + h(s)) + 2∂ 1 F (t, s) ∂ 1 G(t, s, h(t) + h(s)) + h ′ (t) • ∂ 3 G(t, s, h(t) + h(s)) + F (t, s) ∂ 11 G(t, s, h(t) + h(s)) + 2h ′ (t) • ∂ 13 G(t, s, h(t) + h(s)) . + h ′′ (t) • ∂ 3 G(t, s, h(t) + h(s)) + (h ′ (t)) 2 • ∂ 33 G(t

, s, h(t) + h(s))) .

For the derivatives of the function G, we deduce the following estimates |r n (i, j)| ≤ Cn -h(i/n)-h(j/n) i h(i/n)+h(j/n)-2 + j h(i/n)+h(j/n)-2 +|i -j| h(i/n)+h(j/n)-2 (3.24) ≤ Cn -2h min i 2h min -2 + j 2h min -2 + |i -j| 2h min -2 , i, j ≥ 1, |i -j| > 2.

|∂ 1 G(t,
When |i -j| ≤ 2 we deduce from [BS10, Lemma 1 p.13] that Observing the decomposition

S n - 1 2 2h min = V (B h ; 1) n -E[V (B h ; 1) n ] V (B h ; 2) n -U h n V (B h ; 2) n -E[V (B h ; 2) n ] V (B h ; 2) n + U h n - 1 2 2h min
and in view of Proposition 2.2, it is sufficient to show that var(V (B h ; l) n )

E[V (B h ; k) n ] → 0, k, l = 1, 2 (3.27)
to prove Theorem 2.3. We assume again without loss of generality that q = 0, m = 1 and M h = h -1 {h min } = {x}. Using the notations from the previous subsection together with the inequalities (3.16) and (3.17), we deduce the following lower bound, for n large enough and for ǫ small enough:

E[V (B h ; k) n ] ≥ Γ (1)
n,k (γ, ǫ) + Γ

(2) n,k (γ, ǫ) ≥ C ǫ n 1-2h min (ln n) 1/p . Thus, in view of (3.26) we readily deduce the convergence at (3.27) for any h min ∈ (0, 1), which completes the proof of Theorem 2.3.

l

  (x) (resp. h (p) r (x)) the pth left (resp. right) derivative of h at point x.

  u n,k,p (a) := 1 p βn(a) αn(a) y 1/p-1 e -y dy, v n,k,p (a) := 1 p µn(a) τn(a)y 1/p-1 e -y dy, u ′ n,k,p (b) := 1 p β ′ n (b) α ′ n (b)

≤

  |r n (i, j)| ≤ var B h Cn -2h min .(3.25)We recall the identity cov(Z 2 1 , Z 2 2 ) = 2cov(Z 1 , Z 2 ) 2 for a Gaussian vector (Z 1 , Z 2 ). By (3.24) and (3.25) we immediately conclude thatvar(V (B h ; k) n ) ≤ Cn -4h min +1 n i=1 i 4h min -4 ≤ C        n -4h min +1 h min ∈ (0, 3/4) ln n • n -2 h min = 3/4 n -2h min ∈ (3/4, 1)(3.26)

  s, h(t) + h(s))| ≤ C t h(t)+h(s)-1 + |t -s| h(t)+h(s)-1 ,|∂ 3 G(t, s, h(t) + h(s))| ≤ Cln t • t h(t)+h(s)ln s • s h(t)+h(s)ln |t -s| • |t -s| h(t)+h(s) |∂ 11 G(t, s, h(t) + h(s))| ≤ C t h(t)+h(s)-2 + |t -s| h(t)+h(s)-2 |∂ 13 G(t, s, h(t) + h(s))| ≤ C (1ln t) t h(t)+h(s)-1 + (1ln |t -s|)|t -s| h(t)+h(s)-1 |∂ 33 G(t, s, h(t) + h(s))| ≤ C ln 2 t • t h(t)+h(s) + ln 2 s • s h(t)+h(s) + ln 2 |t -s| • |t -s| h(t)+h(s) , which hold for t, s ∈ (0, 1] with t = s and the third inequality holds whenever h(t) + h(s) = 1 (if h(t) + h(s) = 0 we simply have ∂ 11 G(t, s, h(t) + h(s)) = 0). Similar formulas and bounds are obtained for other second order derivatives of R h . Using the boundedness of functions F , h and its derivatives, together with the above estimates and (3.23) we obtain the inequality

See[START_REF] Stoev | How rich is the class of multifractional Brownian motions? Stochastic Processes and their Applications[END_REF] and [ST94, Chapter 6] for more details on Gaussian complex measures.


F Millour
email: fmillour@oca.eu

P Berio

M Heininger

K.-H Hofmann

D Schertl

G Weigelt

F Guitton

W Jaffe

U Beckmann

R Petrov

F Allouche

S Robbe-Dubois

S Lagarde

A Soulain

A Meilland

A Matter

P Cruzalèbes

B Lopez

Data reduction for the MATISSE instrument

Keywords: MATISSE, Interferometry, Data Reduction Software

We present in this paper the general formalism and data processing steps used in the MATISSE data reduction software, as it has been developed by the MATISSE consortium. The MATISSE instrument is the mid-infrared new generation interferometric instrument of the Very Large Telescope Interferometer (VLTI). It is a 2-in-1 instrument with 2 cryostats and 2 detectors: one 2k × 2k Rockwell Hawaii 2RG detector for L&M-bands, and one 1k × 1k Raytheon Aquarius detector for N-band, both read at high framerates, up to 30 frames per second. MATISSE is undergoing its first tests in laboratory today.

INTRODUCTION

MATISSE 1 produces 4-telescopes interferences which are dispersed onto an infrared detector with a spectrograph, simultaneously in the L&M-bands (3.5-4.1µm & 4.6-5.5µm, respectively), and in the N-band (8-13µm). Similarly as AMBER, [START_REF] Millour | Data reduction for the AMBER instrument[END_REF] MATISSE uses an all-in-one multiaxial combination with (SI-PHOT mode) or without (HIGH-SENS mode) photometric separation (see Fig. 1), but due to the limits of the P2VM algorithm that were found for the AMBER instrument, we chose to use the classical Fourier processing instead of a P2VM algorithm for the data reduction software (DRS).

The MATISSE interferogram i(x, λ, t) is governed by the interferometric equation, describing the signal: i(x, λ, t) =

4 i=1 P i + 4 i=1 4 j=i+1 2V ij • P i • P j • cos (2π • f ij • x + Φ ij ) + B (1)
Most, if not all, terms of this equation depend on the coordinates x (position on the detector in the space direction) and λ (position on the detector in the spectral direction), and also time t. However, for clarity, we exhibit their dependence only in each terms description: P i (x, λ, t) is the flux received on the detector from the telescope i. V ij (λ, t) denotes the visibility of the observed object on the telescope pair ij (constant but wavelengthdependent) multiplied by an instrumental and atmospheric transfer function (an unknown, but supposedly slowly variable gain ≤ 1), and Φ ij (λ, t) is composed of the phase of the object (constant but wavelength-dependent) plus the atmospheric phase errors (depending both on λ and t). Let b ij be the combining baseline, i.e. the separation of each of the 6 telescopes pair as seen by the detector, fixed by construction. The fringe rate, or frequency, of the fringe pattern is directly linked to it by f ij (λ) = bij λ . In the case of MATISSE, we have b ij equal to 3 D, 6 D, 9 D, 12 D, 15 D and 18 D, D being the dimension of the output pupil of the instrument, right in front of the detector. Finally, B(x, λ, t) is a dominant and highly variable background contribution from the warm and turbulent atmosphere, the warm optical train of the VLTI and of MATISSE itself. B(x, λ, t) can be up to 10 5 times larger than the flux from the star, and we will see that mitigating this background contribution is a serious task in the MATISSE concept and DRS (as it is for any mid-infrared instrument).

THE DATA REDUCTION SOFTWARE

The MATISSE DRS (called drsmat) aims at producing science-grade spectra, visibilities and phases out of the MATISSE raw datasets. The description of the main algorithms were described in details in several technical documents, and we report here the selected data reduction strategy and algorithms that were implemented in drsmat.

It is implemented in ANSI C and is fully integrated into the ESO pipeline software (i.e. it is interfaced with EsoRex 3 and Reflex 4 ). The implementation uses the ESO-developped Common Pipeline Library (CPL), [START_REF] Mckay | The common pipeline library: standardizing pipeline processing[END_REF] the FITS format for intermediate frames, and the OIFITS format [START_REF] Pauls | A Data Exchange Standard for Optical (Visible/IR) Interferometry[END_REF][START_REF] Duvert | OIFITS 2: the 2nd version of the Data Exchange Standard for Optical (Visible/IR) Interferometry[END_REF] for the final data products (the MATISSE DRS will fully support the OIFITS version 2 format). The Reflex interfaces are made with xml files and also Organisation, Classification, and Association rules (OCA rules), and the GUI part is coded in python. The consortium actively supports ESO in the integration of drsmat in a Reflex workflow that will be used both at ESO and at the consortium sides for the reduction of the MATISSE data.

Overview

The basic input data are the raw frames of the MATISSE instrument, stored as FITS files. Let be i raw (X, Y, t) this raw data intensity at time t, given at the (X, Y ) pixel detector index. Compared to the interferogram i(x, λ, t), the raw data is affected by the following effects:

• Detector bias,

• Detector bad pixels,

• Detector + instrument flat field, plus eventually non-linear effects,

• Instrument distortion map, including a non-linear spectral dispersion law, All these effects must be removed from the data in order to get the interferogram of eq. 1, out of which we extract the MATISSE observables. In addition, a calibration process leads to science-grade data. To do so, the software is organized into a data reduction cascade (shown in Fig. 2) that enables one to produce the calibration maps (mat cal det, mat est flat, mat est shift, mat est kappa), reduce the raw data into oifits files (mat raw estimates), calibrate the data (mat cal oifits), and then reconstruct an image out of the reduced data (mat cal imarec).

The very first step is to estimate the calibration maps like flat field map, bad pixels map, or distortion map out of calibration frames obtained during daytime using the internal calibration source of the instrument. They are used to correct the basic cosmetic of the detector and some optical effects of the instrument.

Detector-specific cosmetics 2.2.1 Bad Pixels, Flat Field

The bad pixel map is derived from a series of cold dark raw frames obtained by increasing the integration time (DIT) from the MINDIT to five times the MINDIT with a minimum of 50 frames per exposure:

• First, spurious cosmics are detected by analyzing the intensity values of each pixel,

• Then, a robust median, a mean, and a standard deviation are computed for each pixel.

• For quality control purposes, the median, mean, and standard deviation of each pixel median of each series are calculated.

• For each pixel, a straight line is fitted to the 5 different median values. This fit results in an offset (pixel bias), a slope (dark current), and a fit quality.

• A pixel is marked as a bad pixel if the fit quality or the slope exceeds a certain limit (N sigmas of the standard deviation of the fit quality or slope for all detector pixels).

The flatfield map is derived from a series of illuminated frames ranging from detector noise up to detector saturation:

• First, the median of the matching series of cold dark frames is subtracted.

• Each frame is processed like the series of cold dark frames, i.e., mean, median, and RMS are computed.

• After normalization, the slopes of the straight lines from the flatfield images are used to compute the master gain map. However, in order to get a "detector" flatfield, the non-flat illumination must be taken into account. Therefore, the implemented algorithm must be adapted to the instrument optics.

The non-linearity map and the flat field are derived by fitting higher-order polynomials to the exposed frames, and saved as maps into fits files. We differentiate two types of flat fields in MATISSE: the detector flat field, which is obtained by a special device put into the beam, which diffuses the light onto the detector in a roughly uniform way, and the observing flat field, which is obtained through the spectrograph with a white lamp or a star without spectral features. The first one enables us to correct for detector pixel-to-pixel gain variations, while the second one allows us to correct for optics transmission variations.

Excess low frequency noise

During the Aquarius detector tests, and thanks to the ESO experience with VISIR, a strong temporal noise was effectively detected in the data. A set of exposures with increasing DIT was taken. This so called Excess Low Frequency Noise (ELFN) showed up in the temporal power spectra, as presented in Fig. 3. Pure white noise would result in a nearly flat temporal power spectrum. The increase at lower frequencies shows the ELFN.

Further investigations have lead to the following: the Aquarius ELFN is not a simple 1/f noise, the origin of the ELFN seems to be inside the semiconductor diode (pixel), the ELFN of different pixels is not correlated, the ELFN could not be reduced by using different readout modes, different multiplexer clocking schemes and delays or different bias voltages, the ELFN is proportional to the detector output signal level, ELFN cannot be seen in dark images, therefore CDS does not help, and finally ELFN was found on other MIR detectors (Si:As, InSb, etc.).

We expect this effect to be negligible for the coherent flux estimate (eq. 11), but likely present in the photometric estimates of MATISSE (eq. 7). The mitigation of this effect will need a fast chopping frequency like in VISIR, whose speed need to be addressed on sky.

Optics calibration 2.3.1 Camera & spectrograph image distortion

Once the bad pixel map and flat field map have been produced, one can continue the instrument calibration with the determination of the image distortion. Since MATISSE contains a long-slit spectrograph, the distortion in the wavelength direction translates into a spectral dispersion law which is non-linear. We chose to treat the problem globally by computing a "shift map".

The shift map is derived from a series of frames containing either a spatial grid (3 separate holes in the slit direction) or a spectral grid introduced by carefully-chosen plastic foils (for the wavelength direction), as shown in Fig. 4. This distortion estimate will be done at the same frequency as the flat field map or bad pixel map estimates.

The distortion can be expressed as the relation between the pixel coordinates (X, Y) and the spectrograph coordinates (x, λ) on the detector: (X, Y ) are the pixel and (x, λ) are the position and spectral dispersion, respectively. This can be written as a matrix relation:

X = a(x, λ)x + b(x, λ)λ Y = c(x, λ)x + d(x, λ)λ (2)
X Y = a(x, λ) b(x, λ) c(x, λ) d(x, λ) x λ = (A) (x, λ) x λ (3)
What we need then is to invert the matrix (A)(x, λ) = a(x, λ) b(x, λ) c(x, λ) d(x, λ) to get the distortion relation:

x λ = (A) -1 (X, Y ) X Y (4)

Spatial Direction:

The spatial grid is first used. The a parameter is estimated by fitting a polynomial to the position of the different features detected on the beams (3 features per beam):

a(x, λ) = a 0 (λ) + a 1 (λ) × x + a 2 (λ) × x 2 + ... (5)

Spectral direction:

We are then only interested in determining the λ-dependent distortion. Therefore, a "normal" λ-calibration frame (i.e. with the set of absorption plastic foils, or using sky absorption lines) should be available for that. In that case, one can determine the d coefficient:

d(x, λ) = d 0 (x) + a 1 (x) × λ + d 2 (x) × λ 2 + ... (6)
We note here that this first step will give a coarse wavelength-calibration, which will be refined in a further step using telluric lines.

Shift and zoom, κ-matrix

As shown in Fig. 1 and Fig. 4, the MATISSE photometric beams have a different width than the interferometric beam, for obvious pixels-saving reasons. Therefore, since the spatial filtering is not perfect (due to the fact that we use pinholes instead of fibers, for technology readiness reasons), one needs to scale the photometric beams to the interferometric beam size by shifting and zooming them to the right place and scale.

In addition, one has to get the flux ratio coefficient between the flux measured in the interferometric beam and the one measured in the photometric beam, in order to have all the photometric information available to compute the observables.

The Shift and Zoom coefficients, together with the κ-matrix are computed using 4 illuminated frames, one for each of the telescope beams, the other beams being closed by shutters (see Fig. 5). The source is preferably artificial, part of MATISSE instrument, or could be a bright unresolved astronomical target as well. photometric contribution P ′ i (x, λ) of the telescope i to the interferometric beam can be expressed as following:

P ′ i (x, λ, t) = κ ij [x ′ , λ, s j (λ), z j (λ)] × P ij (x ′ , λ, t) (7)
where P ij (x ′ , λ, t) is a vector containing the contribution of the telescope i to the photometric beams j (after chopping correction, see section 2.4.1), κ ij is the linear transformation matrix of the intensities in the photometric channels into the interferometric channel (the so-called "κ-matrix"), and s j (λ), z j (λ) are the shift offset and zoom coefficient to match the photometric beam j into the interferometric beam.

We determine therefore the κ-matrix and the shift-and-zoom coefficients by fitting the photometric beam shape and intensity to the interferometric beam for each selected wavelength.

Applying cosmetics

The first steps of the data reduction are:

How to tackle the multiple-stage modulations in MATISSE

As described in Petrov et al. (2016), 8 MATISSE includes a three-stage modulation process for the mitigation of the sky background (as MATISSE works in the mid-infrared). These stages are namely chopping, spatial modulation and temporal modulation, in addition to beam commutation, which is used to clean up the phases out of the instrumental effects.

Chopping

Chopping produces alternatively science sky frames with a switching frequency as fast as 1Hz. The first step is to remove these sky frames to the science frames, resulting in removing the bulk of the background contribution (≈ 99%):

i(x, λ, t) = i obj (x, λ, t) -i sky (x, λ, t) (8)

Spatial modulation

Spatial modulation separates, in the Fourier plane, the fringes energy ("fringe peak") from the flux energy ("photometric peak"). Therefore, a simple Fourier processing along the x axis will cancel further out the contribution of the background to the data. Correlated fluxes are computed just by Fourier-Transforming the clean interferogram:

I ′ (u, λ, t) = FT x [i(x, λ, t)] = F 00 (u, λ, t) + 4 i=1 4 j=i+1 F ij (u -f ij , λ, t) • V ij (λ) • e ıΦij (λ,t) (9)
where u is a spatial frequency, f ij is the fringe rate, as explained below eq. 1, V ij (λ) is the target visibility, Φ ij (λ, t) is the phase of the target plus piston, and F ij (u, λ, t) is a function of the photometric contributions P ′ i (x, λ), involving Fourier Transforms and correlations, representing the shape of the peaks. We note that u F ij (u, λ, t) =

x P i (x, λ, t) • P j (x, λ, t), and that F 00 (u, λ, t) represents the photometry peak.

The background residuals left after chopping contribute only to F 00 (u, λ, t) at first order, as the background does not contain any coherent information. At the frequency of each fringe peak f ij = b ij /λ, we have:

I ′ ij (λ, t) = I ′ (f ij , λ, t) = F ij (0, λ, t) • V ij (λ) • e ıΦij (λ,t) + B res ij (λ, t) (10)
where the contribution of the residual background B res ij (λ, t) = F 0 (f ij , λ, t) is much smaller than at the frequency 0 (as seen in Fig. 7 on top and bottom row).

Temporal modulation

In addition to the spatial modulation, a temporal modulation is introduced by piezoelectric actuators, changing the optical path difference (OPD) ∆ mod ij (t) by a fraction of the wavelength between consecutive frames. Let us consider that there are N T steps over an OPD range of ∆ ij during one modulation cycle T . The modulation is then ∆ mod ij (t) = k NT ∆ ij , k being the modulation step number. On needs to consider the fringes "freezed" over the modulation time T in order to consider ∆ atm ij (t) constant. Therefore, the modulation cycle must be done over one atmosphere coherence time T ≈ τ 0 , or in conjunction with a fringe tracker which freezes the fringes over a time T ≈ τ FT .

A de-modulation is applied to the previous Fourier-Transform, baseline-by-baseline, and we integrate the fringe peak I ′ ij (λ, t) over one modulation cycle (the N T steps of modulation). We do so by multiplying it with a phasor containing the counter-modulation:

I ′′ ij (λ, T ) = 1 N T NT k=0 I ′ ij (λ, t) × e -2iπ∆ mod ij (t)/λ = F ij (0, λ, t) • V ij (λ) • e iΦij (λ,t) + δB res ij (λ, t) (11)
The background residual is further reduced by the modulation, leading to the final residual background δB res ij (λ, t):

res ij (λ, t) = B res ij (λ, t) × 1 N T NT k=0 e [-2iπk•∆ij /(NT •λ)] (12)
The effect of demodulation can be seen in Fig. 7, middle row. The combination of spatial modulation, temporal modulation and chopping are the three ways selected in MATISSE to tackle the dominance of the background. Please note that neither the spatial nor the temporal modulation can be applied to the photometries. Therefore, the rejection of the sky background in the photometric estimates rely only on Chopping, and will be of much lower quality.

Estimators

Once the complex coherent flux I ′′ ij (λ, T ) has been computed, it is possible to compute the estimators, i.e. the raw visibilities and phases of the object. These estimators are split between incoherent estimators (speckle-like), that can be computed without the knowledge of the optical path difference, and coherent estimators that need an appropriate atmospheric OPD estimate, be it simply due to the turbulence or to higher-order effects (longitudinal dispersion).

"Incoherent" estimators

Squared visibility Squared visibilities are computed as follows:

V 2 ij (λ) = u ′ (u, λ, t)| 2 -β t 2 • x P a ′ (x, λ, t) • P b ′ (x, λ, t) t ( 13
)
where u is the frequency around the considered baseline, and β a bias estimated on I ′ outside the range of frequencies where the fringe peaks are present, and P a ′ (x, λ, t) and P b ′ (x, λ, t) are the estimates of the photometric fluxes transformed into the interferometric channel (see eq. 7).

In the case where there are no photometries recorded (like in HIGH-SENS mode without a photometric acquisition sequence), the squared visibilities are simply not computed. Instead, the user should use the correlated fluxes, i.e. the upper part of eq. 13:

CF 2 ij (λ) = u |I ′ (u, λ, t)| 2 -β t ( 14
)
Closure Phase The bispectrum is calculated by:

O (3) ij (λ) = I ′ (u, λ, t) × I ′ (v, λ, t) × I ′ * (u + v, λ, t) t ( 15
)
where u, v, and u+v are the frequencies of the baselines involved in the closure relation. An example MATISSE bispectrum can be seen in Fig. 8. Out of the 5 peaks that can be seen, just 4 contain useful information from the 4 closure phases available at 4 telescopes (the middle-left peak is a spurious peak coming from the contamination of the central photometric peak). This average bispectrum contains an additive noise bias terms plus a multiplicative one, because of the all-inone combination of the interferograms. These biases need to be subtracted in order to get closure phases without systematic errors, following specific recipes. 9, 10

Figure 1 .

 1 Figure 1. Expected (Left) vs. real data (Right) beams implementation on the MATISSE Hawaii detector (L-band). The wide central stripe contains the dispersed fringes (vertical spectral dispersion), while the 4 narrow stripes on the side contain the 4 single-telescope photometry used for flux calibration.

Figure 2 .

 2 Figure 2. The MATISSE drsmat pipeline cascade, showing all the steps and files necessary to obtain a final science-grade dataset (calibrated oifits and/or reconstructed images).

Figure 3 .

 3 Figure 3. Average temporal power spectrum at several exposure times and increasing illumination. The ELFN shows up at frequencies lower than 1 Hz.

Figure 4 .

 4 Figure 4. Three panels on the left: Spatial grid, made with three pinholes at different positions. Note that the photometric beams position is inverted compared to the interferometric beam one. Right: Spectral grid, made with plastic foils with spectral features. These snapshots were recorded in the N-band.

Figure 5 .

 5 Figure 5. Kappa matrix acquisition sequence in the N-band, with alternatively just one out of four open shutter.

1 .

 1 Subtracting the average (cold) dark DARK(X, Y ) in each frame, 2. Compensating of the space-variant gain in each frame by division of each frame through the FFM(X, Y ) (flat-field map), 3. Interpolating detector bad pixels in each frame with the BPM(X, Y ) (bad pixel map), 4. Applying the distortion matrix to transform (X, Y ) to (x, λ), the coordinates of the data, 5. Adjusting the photometric contribution in the interferometric beam by applying the κ-matrix, shift, and zoom coefficients The result of treating the three first effects is shown in Fig. 6. The two last corrections are being tested right now on the mounted instrument in the Nice laboratory. The results of these processes are: a cleaned-up fringe pattern i ′ (x, λ, t), and cleaned-up photometric estimates P ′ i (x, λ, t).

Figure 6 .

 6 Figure 6. Left: The raw image from the Aquarius detector (7-13 µm) contains four photometric channels, the interferometric channel and 16 reference areas, showing bad pixels and detector channel offsets. Right: The same image after the calibration step contains the photometric channels and the interferometric channel only.

Figure 7 .

 7 Figure 7. Left column: Fourier transform of the MATISSE fringe pattern. Right column: cut at one wavelength. Top: Simulated fringe pattern after the chopping correction. Middle: Same as before, but after the temporal demodulation applied. Bottom: Same as first one but on the real first fringes data shown in Fig. 1. The (de)modulation is not yet operational on the instrument.

Figure 8 .

 8 Figure 8. Example bispectrum values computed on simulated MATISSE data. Top graph: real part, and bottom graph: imaginary part. The middle-left peak is spurious.

ACKNOWLEDGMENTS

We are grateful to ESO, CNRS/INSU, and the Max-Planck Society for continuous support in the MATISSE project.

"Coherent" estimators

Coherent estimators need an appropriate OPD estimate. We describe here briefly what we envision to include in drsmat. The complex coherent flux is provided by the equation 11.

The phase term Φ ij (λ, t) depends on the variations of the atmospheric properties (temperature, pressure), producing an achromatic variation of OPD (but not phase), and it depends also on the refractive index n(λ)) of air, which is chromatic and dependent of its water vapor composition. It can be expressed as:

the different terms are the following:

, chromatic term with a(λ) wavelength-dependent but roughly static, [START_REF] Ciddor | Refractive index of air: new equations for the visible and near infrared[END_REF][START_REF] Mathar | Refractive index of humid air in the infrared: model fits[END_REF] and b(t) wavelength independent, but with strong variations in time 13

ij /λ+Φ diff ij (λ): the object phase, which can be decomposed in a Taylor series containing the differential phase Φ diff ij (λ).

Achromatic Optical Path Difference estimates: MATISSE will compute the achromatic OPD using a VEGA-like algorithm: 14, 15 a 2D-Fourier Transform (FT) is computed, to take profit of the multi-axial beam combination of MATISSE. The fringe peak in the power spectrum appears as an offset peak, whose offset is proportional to the OPD. The great advantage of this method is that the peak is concentrated onto a small location of the 2D map, whereas the noise is spread everywhere, hence contributing to a smaller fraction of the peak noise. In addition, there is no ambiguity at OPD 0, contrary to a MIDI-like estimate.

This 2D-FT can be computed in two steps, first in the spatial direction and second in the spectral direction, allowing us to counteract the OPD modulation pattern before elevating the Fourier Transform to its square, and hence increasing tremendously the SNR of the peak detection. This is also the way the fringe pattern will be detected in real time with the instrument.

Another option implemented for OPD determination is a fit to the complex coherent flux phase, in the very same way as it is done for AMBER. [START_REF] Tatulli | Interferometric data reduction with AMBER/VLTI. Principle, estimators, and illustration[END_REF] The computed atmospheric OPD contains then the two terms: OPD = ∆ atm ij (t) + Φ

(1) ij .

Chromatic OPD estimates: ESO provides information through the ambient conditions monitor that will be sufficient to correct the data from the chromatic dry and wet air dispersion, down to accuracies of 3 degrees in K band and down to 10 degrees in N-band.

The static water vapor term a(λ) will be computed using the pressure P , temperature T (in the tunnels), and relative humidity of the ambient air R h . We add to these the partial pressure of CO 2 , which can be hardcoded to an average value at Paranal. The fraction of humidity is converted into partial pressure of water vapor p wv = R h × P . The λ-dependent shape of the chromatic phase a(λ) is then be determined based on computed optical index of air n as a function of wavelength. [START_REF] Ciddor | Refractive index of air: new equations for the visible and near infrared[END_REF][START_REF] Mathar | Refractive index of humid air in the infrared: model fits[END_REF] The next step is to compute phases, and therefore we also use the information of the delay lines position to compute the difference of delay for each baseline: if A i is the static optical path introduced by the delay line i, and OP L i the dynamic path, A i + OP L i is the total path of air introduced by each delay line. Per baseline ij, the path to use for the chromatic OPD is therefore

The chromatic phase is then computed by φ chrom (λ) = 2πa(λ) δij λ , with n the computed index of air, and λ the wavelength of interest.

The water vapor phase offset, time-variable, ∆ atm ij (t) • b(t), will be obtained by averaging the phase term in the wavelength-direction after subtraction of the atmospheric OPD. As a consequence, it also removes the object-related phase offset Φ

Differential phase: After that step, the phase is transformed into a differential phase the standard way. [START_REF] Millour | Differential Interferometry with the AMBER/VLTI instrument: Description, performances and illustration[END_REF][START_REF] Millour | Circumstellar Matter Studied by Spectrally-Resolved Interferometry[END_REF][START_REF] Millour | Interferometry concepts[END_REF] We recall it briefly here. First, we multiply the Fourier transform of the correlated interferograms by the counter OPD phasor. This yields the OPD-free complex coherent flux:

Averaging this OPD-free complex coherent flux over time yields the following products:

: correlated flux degraded by the not corrected beam overlap

• Φ diff ij (λ) : differential phase, with any wavelength-linear phase term set to zero.

Coherent visibility: Finally, it is possible to continue the processing to lead to a linear visibility estimate.

For that, we divide the correlated flux by the photometric factor. This is the last optional step, yielding the linear visibility estimate:

In the case where there are no photometries recorded (like in HIGH-SENS mode), this term is not computed, as for the squared visibilities.

Data calibration

To calibrate the data, we apply the same recipes as in AMBER, [START_REF] Millour | Advanced" data reduction for the AMBER instrument[END_REF] i.e.:

• We start with searching calibrator stars diameters in published catalogs,

• we divide the calibrators visibilities by their expected visibilities and store the transfer function result in new files,

• we interpolate the transfer function to the time of the science observations with a Gaussian-weighted average (typical FWHM 1hr),

• we divide the science visibilities with the interpolated transfer function and store the result as calibrated science-grade data files.

drsmat uses the ESO-developed Reflex interface to run all the steps, and the calibration part contains a display of the result, allowing the user to select the relevant calibration stars (see Fig. 9)

Image reconstruction

The last step of drsmat is the image reconstruction. MATISSE will use the IRBIS algorithm, which was specially developed for it. The content of IRBIS is extensively described in another paper of this conference. [START_REF] Hofmann | Image reconstruction method irbis for optical/infrared long-baseline interferometry[END_REF] The drsmat software displays the results of the reconstruction through a python interface, shown in Fig. 9.

CONCLUSION

We presented the MATISSE data reduction software drsmat as it has been implemented for the MATISSE instrument. drsmat produces calibrated dispersed visibilities, closure phases and differential phases from the beginning, and also includes the IRBIS algorithm for image reconstruction. MATISSE is presently undergoing laboratory tests to verify that it meets all the necessary requirements to reach its science goals. The tests of drsmat is a significant part of these tests.


F Millour
email: fmillour@oca.eu

J Hron

A Chiavassa

G Weigelt

A Soulain

Z Khorrami

A Meilland

N Nardetto

C Paladini

A Domiciano De Souza

G Niccolini

K.-H Hofmann

D Schertl

P Stee

P Bendjoya

F Thévenin

F Vakili

P Berio

T Lanz

A Matter

P Cruzalèbes

R Petrov

B Lopez

Grown-up stars physics with MATISSE

MATISSE represents a great opportunity to image the environment around massive and evolved stars. This will allow one to put constraints on the circumstellar structure, on the mass ejection of dust and its reorganization, and on the dust-nature and formation processes. MATISSE measurements will often be pivotal for the understanding of large multiwavelength datasets on the same targets collected through many high-angular resolution facilities at ESO like sub-millimeter interferometry (ALMA), near-infrared adaptive optics (NACO, SPHERE), interferometry (PIONIER, GRAVITY), spectroscopy (CRIRES), and mid-infrared imaging (VISIR). Among main sequence and evolved stars, several cases of interest have been identified that we describe in this paper.

HOW MANY EVOLVED STARS CAN BE OBSERVED WITH MATISSE?

First of all, we need to emphasize that MATISSE will mainly focus on dusty stars, as the mid-infrared is the perfect match to collect data on the dust quantity and composition around a given object.

To illustrate that, we selected a few topics of interest that will be developed in this paper, and counted the number of stars that will be observable with VLTI/MATISSE. We included a quantity of non-dusty targets (regular WR stars and Be stars) in the sample to compare the performances of MATISSE with its prime targets.

To do so, one need to take into account the sensitivity of the instrument itself in the L-and N-bands, [START_REF] Matter | Matisse specifications[END_REF] of course, but also the VLTI infrastructure subsystems sensitivity, which can be sometimes the limiting factor, especially for red targets. We need therefore to take into account: the telescope guiding limit in the V-band (V=13.5 and V=17 for ATs STRAP tip-tilt and UTs MACAO AO, respectively), and the fringe tracker limit in the K-band for GRA4MAT and other devices like IRIS (K=7.5 for ATs and K=10 for UTs).

We considered star lists on Wolf-Rayet stars (WR), Red supergiant stars (RSG), Asymptotic Giant Branch stars (AGB), Be stars and B[e] supergiant stars visible from Paranal ; we extracted them from the following catalogs:

• Rosslowe & Crowther (2014) for WR stars, [START_REF] Rosslowe | Spatial distribution of Galactic Wolf-Rayet stars and implications for the global population[END_REF] • Hoffleit (1991) for RSG stars, [START_REF] Hoffleit | The Bright star catalogue[END_REF] • we scanned through ADS publications on the topic for AGB stars,

• Frémat et al. (2005) [START_REF] Frémat | Effects of gravitational darkening on the determination of fundamental parameters in fast-rotating B-type stars[END_REF] and Yudin (2001) [START_REF] Yudin | Statistical analysis of intrinsic polarization, IR excess and projected rotational velocity distributions of classical Be stars[END_REF] for Be stars,

• For B[e] supergiant stars we used our own catalog of stars.

When available, we retrieved the SED of each object (V, K, L and N magnitudes) and ISO/IRAS spectra (if it exist). Magnitudes of these targets were compared with the theoretical limits of MATISSE, considering spectrally-usable data, i.e. data with a SNR≥ 3 per spectral channel during the expected optimal DIT and exposure time combination, with and without fringe tracker. [START_REF] Matter | Matisse specifications[END_REF] The result of this study is shown in Fig. 1. One can see that, out of our selected catalogs, roughly one third of the considered AGB targets, two third of the supergiant B[e] stars, and half of the considered red supergiant stars can be observed with MATISSE in its low spectral resolution mode (R≈35). We are likely limited by our sample size in these cases.

On the other hand, very red targets like dusty Wolf-Rayet stars, or very blue targets like naked Wolf-Rayet stars or Be stars, will be more difficult to observe with MATISSE due to either the sensitivity limits of VLTI in the V band (very red targets), or to the MATISSE sensitivity limit (very blue targets). In both cases, the adjunction of a K-band fringe tracker to MATISSE, like the foreseen GRA4MAT project of ESO, will improve the situation by a large factor (especially when using the ATs).

That matter of fact is even more striking when using the high-spectral resolution, where we will be able to resolve the kinematic motion of gaz in the circumstellar envelopes. From a handful of targets that will be observable with MATISSE alone, the adjunction of an external fringe tracker will allow the instrument to go at full throttle for the delivery of scientific results.

MASSIVE STARS

Massive stars are stars with a mass above 8 solar masses (M ). They evolve fast and blast as supernovae at the end of their lives, not without having gone through a fireworks of different physical processes. Among them, one can count the main sequence O and B-type stars, red, yellow and blue supergiant stars, Wolf-Rayet stars, and other types of exotic stars. We invite the reader to take a look to the different other papers of this conference for more details. 6, 7

Dusty blue supergiant stars

Several blue supergiant stars (Wolf-Rayet, Luminous Blue Variable, supergiant B[e] stars) produce large amounts of dust before exploding as supernovae. Though they are not the main dust-producers in our Galaxy, they do vastly influence their surroundings by inputting kinetic energy and dust in their vicinity. In some cases, they can even trigger star formation.

The circumstellar environment (CE) around blue supergiant stars is often dense and highly anisotropic. This anisotropy could be driven either by fast-rotation and/or by the presence of a companion star. For example, fast-rotating B and O stars will present a denser and slower wind in the equatorial direction, and a lighter and faster wind at the poles, according to the usual sketch. [START_REF] Lamers | The formation of outflowing disks around early-type stars by bi-stable radiation-driven winds[END_REF] At the equator, the density can be high-enough to veil the UV radiation from the star, allowing dust to condense into a disk-like structure, visible in the infrared. Such a break of symmetry of the CE has strong implications on the later phases of stellar evolution, as it can affect the stellar angular momentum and act as a brake for ejected material.

On the other hand, a secondary star orbiting around an O or a B star can affect the wind of the main star through gravitational influence, and concentrate material in the orbital plane. The consequence is similar as for rotating stars, i.e. an over-density in the equatorial plane that can trigger dust formation. However, the shape of such a disk-like structure would likely be different than in the first case, and there are opportunities of detecting the companion star, depending on the flux ratio between the two stars. [START_REF] Millour | A binary engine fuelling HD 87643's complex circumstellar environment. Determined using AMBER/VLTI imaging[END_REF][START_REF] Kraus | On the Nature of the Herbig B[e] Star Binary System V921 Scorpii: Geometry and Kinematics of the Circumprimary Disk on Sub-AU Scales[END_REF] To prepare the MATISSE program on supergiant B[e] stars, we simulated what are the imaging performances of MATISSE in the case of a disk-like structure and a possible influence of fast-rotation of the hot star. To prepare and interpret the MATISSE observations we used a NLTE Monte Carlo radiation transfer code to compute a grid of supergiant B[e] models dedicated to mid-and near-IR interferometry. [START_REF] Domiciano De Souza | Grid of Supergiant B[e] Models from HDUST Radiative Transfer[END_REF] The case of binarity is expected to be much easier, as it was already demonstrated that imaging with few baselines can be done in such a case. [START_REF] Millour | A binary engine fuelling HD 87643's complex circumstellar environment. Determined using AMBER/VLTI imaging[END_REF][START_REF] Kraus | On the Nature of the Herbig B[e] Star Binary System V921 Scorpii: Geometry and Kinematics of the Circumprimary Disk on Sub-AU Scales[END_REF] We simulated MATISSE with homemade scripts and produced OIFITS 12 files. We input these files into the MiRA image reconstruction software developed by E. Thiébaut. [START_REF] Thiebaut | Image reconstruction in optical interferometry[END_REF] The results of our study is shown in Fig. 2. One can see that we are able to reconstruct most of the features seen in the model images: the central point source, the extent of the disk, and the inner hole seen in the 60 degrees inclination images. This study confirms that MATISSE will be a perfect imaging machine for these types of stars.

The case of η Carinae

The Luminous Blue Variable (LBV) phase is a short-lived phase in the life of massive stars. The strong mass loss in this phase is not well understood. The LBV η Car provides a unique laboratory for studying the evolution of massive stars and the massive stellar wind during the LBV phase. η Car is also an eccentric binary with a period of 5.54 yr. X-ray studies demonstrate that there is a wind-wind collision zone, which is depends on the orbital phase. Infrared interferometry with the VLTI allows us to study the stellar wind of the primary star and the wind-wind collision zone of the primary and secondary star with unprecedented angular resolution of about 5 mas and simultaneously with high spectral resolution. Using the VLTI, the diameter of η Car's wind region was measured to be about 4.2 mas in the K-band continuum [START_REF] Van Boekel | Direct measurement of the size and shape of the presentday stellar wind of eta Carinae[END_REF][START_REF] Kervella | Interferometric observations of η Carinae with VINCI/VLTI[END_REF][START_REF] Weigelt | Near-infrared interferometry of ¡ASTROBJ¿η Cari-nae¡/ASTROBJ¿ with spectral resolutions of 1 500 and 12 000 using AMBER/VLTI[END_REF] (50% encircled intensity diameter). In the He I 2.059 µm and the Brγ 2.166 µm emission lines, diameters of 6.5 and 9.6 mas were measured, respectively. [START_REF] Weigelt | Near-infrared interferometry of ¡ASTROBJ¿η Cari-nae¡/ASTROBJ¿ with spectral resolutions of 1 500 and 12 000 using AMBER/VLTI[END_REF] For comparison, the radius of the primary star of η Car is on the order of 100 solar radii ∼ 0.47 au ∼ 0.20 mas. [START_REF] Hillier | On the Nature of the Central Source in η Carinae[END_REF] This means that the measured 4.8 mas Brγ wind radius is about 25 times larger than the stellar radius. The measured line visibilities agree with predictions of the radiative transfer model of Hillier et al. [START_REF] Hillier | On the Nature of the Central Source in η Carinae[END_REF] VLTI-AMBER observations (performed at the beginning of 2014, about 5 to 7 months before the August 2014 periastron passage) allowed the reconstruction of velocity-resolved aperture-synthesis images in more than 100 different spectral channels distributed across the Brγ 2.166 µm emission line. [START_REF] Hofmann | Image reconstruction method irbis for optical/infrared long-baseline interferometry[END_REF][START_REF] Weigelt | VLTI-AMBER velocity-resolved aperture-synthesis imaging of eta Carinae with a spectral resolution of 12 000[END_REF] The intensity distribution of the obtained images strongly depends on wavelength. Interestingly, in the blue wing of the Brγ 2.166 µm emission line, the wind region is much more extended than in the red wing. At radial velocities of approximately -140 to -376 km/s measured relative to line center (i.e., in the blue line wing), the intensity distribution is fan-shaped with a position angle of the symmetry axis of ∼126 degree. The fan-shaped structure extends approximately 8 mas (∼19 au) to the south-east and 6 mas (∼14 au) to the north-west (measured at the 16% intensity contour). Three-dimensional smoothed particle hydrodynamic simulations of η Car's colliding winds predict a very large wind density distribution that is more extended than the undisturbed primary star wind density distribution. At the time of our observations, the secondary star of η Car was in front of the primary star, and therefore the extended wind collision cavity was opened up into our line of sight. The comparison of the model images with the shape and wavelength dependence of the intensity distributions of the reconstructed VLTI images suggests that the obtained VLTI images are the first direct images of the innermost wind-wind collision zone.

In the red line wing (positive radial velocities), the extension of the resolved wind structure is much smaller than in the blue wing, because the red-wing light is emitted from the wind region behind the primary star, where the primary wind is not disturbed by the wind collision (i.e., at positive radial velocities, we are looking through the optically thin wind collision zone in front of the primary). The η Car channel maps provide velocitydependent image structures that can be used to test three-dimensional hydrodynamical models of the massive interacting winds.

Dust plumes around Wolf-Rayet stars

We also investigated the potential of MATISSE to detect asymmetries in the dusty environment of massive blue stars, with an application to the detection of several spiral dust plumes around dusty Wolf-Rayet stars. These dust plumes likely trace the dust forming in the wind collision zone of a binary system composed of a WR and an O star. Many dust shells have been detected around Wolf-Rayet stars, but few of them were resolved into spiral plumes. MATISSE has the potential to image most of these dust shells at a few milliseconds of arc resolution, enabling the detection of many more spiral plumes, leading to a confirmation that WC8 and WC9 spectral types are linked to the inner system binarity.

To illustrate the MATISSE potential, we simulated observations of a phenomenological "pinwheel" model with MATISSE using the readily available ASPRO software that has been updated to take into account the MATISSE-specific noises. The obtained OIFITS files were fed into the MATISSE official image reconstruction tool IRBis, which notably won the interferometry beauty contest of this year. [START_REF] Sanchez | The 2016 interferometry beauty contest[END_REF] The result is shown in Fig. 3, evidencing that MATISSE is indeed able to resolve a "pinwheel" nebula in 3 nights of observation.

Model Image

Young massive clusters

Young massive clusters, are gravitationally bound groups of newly formed stars. These objects provide a rare opportunity to study the formation of massive stars and their evolution in the early stages. Among the wellknown young clusters, R136 is located in the heart of 30Doradus (in LMC) and hosts so far the most massive stars known in the Local Universe. These stars are born in the large HII regions, embedded in gas and dust. MATISSE can be used to resolve the crowded core of R136-like clusters with higher angular resolution than conventional telescopes equipped with adaptive optics (e.g. SPHERE), and at longer wavelengths (L-and N-bands), in order to overcome the problem of interstellar extinction.

We simulated R136 using the publicly available N-body gravitational simulation Nbody6 code † . The left part of Figure 4 shows the inner 0.075 × 0.075 pc of the simulated cluster, which correspond to the FoV of MATISSE using the UTs (0.3"). This snapshot of the core from the numerical Nbody6 simulation is done at the age of 2 Myr in the N-band and at the distance of the LMC (50 Kpc). We used the ASPRO tool from JMMC ‡ to simulate the MATISSE observables (visibilities, closure phases and differential phases), and the IRBis software to reconstruct the synthetic image. The result can be seen in the right panel of Figure 4. The core massive stars have separations down to 24 mas. However, we also put a secondary companion in a binary system with separation of 3.5 mas (seen close to the lower-left star in the left part of Figure 4), and this companion star cannot be resolved by MATISSE in the N band. However, we did not perform image reconstruction tests in the L band, and it could be that this companion star is resolved at this wavelength, where MATISSE has a factor 3 higher angular resolution. Further tests need to be done to fully characterize the MATISSE performances in the context of star cluster observations.

Red supergiant stars

Red supergiants (RSGs) are strong contributors to the chemical and dust enrichment of the Galaxy through their mass-loss (2 × 10 -7 M /yr -3 × 10 -4 M /yr, de Beck et al. 2010) and the coming blast as supernova. Mass-loss process still remains poorly understood. RSGs do not experience regular flares or pulsations and the triggering of their mass loss remains a mystery. It may be linked to the magnetic activity, the stellar convection, and the radiation pressure on molecules and/or grains, or all these together. The stellar surface and circumstellar dynamics and constitution of these stars have to be studied in detail to unveil the mass-loss mechanism. MATISSE will open new insight at wavelengths still poorly observed for RSGs. The observations will range from the inner stellar photosphere (L and M band probing molecular transitions and pseudo-continuum) to the outer enveloppe (N band). The temporal variation is crucial to make a step forward into the understanding of the stellar dynamics and the link with the mass-loss mechanism.

LOWER MASS STARS: AGB & PULSATING STARS

AGB stars

Most of the low and intermediate mass stars end their lives on the Asymptotic Giant Branch (AGB) phase. During this phase, pulsation and radiation pressure on dust leads to a phase of strong mass loss, during which gas and dust enriched by the products of the star's nucleosynthesis will be ejected. This mass loss is crucial for the chemical enrichment of the ISM and therefore for the chemical evolution of galaxies. A pivotal aspect of the mass-loss process is its geometry, i.e. the density distribution of the circumstellar envelope of the AGB stars at different scales and different evolutionary phases.

To understand the mass-loss process, it is essential to study the mass-loss from very deep inside the star up to the interface with the ISM. Due to its broad wavelength coverage MATISSE is the unique instrument to study the different dust and molecular species present in the atmospheres and envelopes of AGB stars. The aim is to achieve a complete view from the upper photosphere to the outer envelope layers and to complement Herschel and MIDI surveys (see Paladini et al., this conference).

Envelopes around Cepheid stars

A powerful way of constraining the period-luminosity (PL) relation is to use the Baade-Wesselink (BW) method of distance determination. The basic principle of the BW method is to compare the linear and angular size variation of a pulsating star in order to derive its distance through a simple division. The angular diameter is either derived by interferometry [START_REF] Kervella | Cepheid distances from infrared long-baseline interferometry. I. VINCI/VLTI observations of seven Galactic Cepheids[END_REF] (hereafter IBW for Interferometric Baade-Wesselink method), or by using the InfraRed Surface Brightness (IRSB) relation. [START_REF] Storm | Calibrating the Cepheid period-luminosity relation from the infrared surface brightness technique. I. The p-factor, the Milky Way relations, and a universal K-band relation[END_REF][START_REF] Storm | Calibrating the Cepheid period-luminosity relation from the infrared surface brightness technique. II. The effect of metallicity and the distance to the LMC[END_REF] However, when determining the linear radius variation of the Cepheid by spectroscopy, one has to use a conversion projection factor from radial to pulsation velocity. [START_REF] Nardetto | Self consistent modelling of the projection factor for interferometric distance determination[END_REF][START_REF] Mérand | The projection factor of δ Cephei. A calibration of the Baade-Wesselink method using the CHARA Array[END_REF] In this method, angular and linear diameters have to correspond to the same physical layer in the star to provide a correct estimate of the distance. [START_REF] Nardetto | High-resolution spectroscopy for Cepheids distance determination. II. A period-projection factor relation[END_REF][START_REF] Nardetto | High-resolution spectroscopy for Cepheids distance determination. V. Impact of the cross-correlation method on the p-factor and the γ-velocities[END_REF] Consequently, the circumstellar environment (CSE) around Cepheids should be considered in both versions of the BW method, particularly when deriving the angular diameter curve. The impact of CSE on the period-luminosity relation of Cepheids need also to be established. Envelopes around Cepheids have been discovered by long-baseline interferometry in the K Band with VLTI 27 and CHARA 28 (see Fig. 5). Since then, four Cepheids have been observed in the N band with VISIR and MIDI [START_REF] Kervella | The circumstellar envelopes of the Cepheids l Carinae and RS Puppis. Comparative study in the infrared with Spitzer, VLT/VISIR, and VLTI/MIDI[END_REF][START_REF] Gallenne | Extended envelopes around Galactic Cepheids. IV. T Monocerotis and X Sagittarii from mid-infrared interferometry with VLTI/MIDI[END_REF] and one with NACO. [START_REF] Gallenne | Spatially extended emission around the Cepheid RS Puppis in near-infrared hydrogen lines. Adaptive optics imaging with VLT/NACO[END_REF][START_REF] Gallenne | Thermal infrared properties of classical and type II Cepheids. Diffraction limited 10 µm imaging with VLT/VISIR[END_REF] Some evidences have also been found using high-resolution spectroscopy. [START_REF] Nardetto | High-resolution spectroscopy for Cepheids distance determination. IV. Time series of Hα line profiles[END_REF] The size of the envelope seems to be at least 3 stellar radii and the flux contribution in K band is from 2% to 10% of the continuum, for medium-and long-period Cepheids respectively, while it is around 10% or more in the N band (estimated also from SED derived with VISIR). Recently, Nardetto et al. (2016, in press) found a resolved structure around the prototype classical Cepheid delta Cep in the visible spectral range using visible VEGA/CHARA interferometric data observations. These data are indeed consistent in first approximation with a quasi-hydrostatic model of pulsation surrounded by a static CSE with a size of 8.9 ± 3.0 mas and a relative flux contribution of 0.07 ± 0.01. A model of visible nebula (a background source filling the field of view of the interferometer) with the same relative flux contribution is also consistent with the data at small spatial frequencies. However, in both cases, they find discrepancies in the squared visibilities at high spatial frequencies (maximum 2σ) with two different regimes over the pulsation cycle of the star, φ = 0.0 -0.8 and φ = 0.8 -1.0. One possibility could be that the star is lighting up its environment differently at minimum (φ = 0.0 -0.8) and maximum (φ = 0.8 -1.0) radius. This reverberation effect would then be more important in the visible than in the infrared since the contribution in flux of the CSE (or the background) in the visible band is about 7% compared to 1.5% in the infrared (Merand et al. 2015).

The processes at work in infrared and in the visible regarding the CSE are different. We expect thermal emission in the infrared and scattering in the visible. MATISSE offers an unique opportunity to study the envelopes of Cepheids. Determining the size of these envelope (as a function of the spectral band) and their geometry (as a function of the pulsation phase), will bring insights in the links between pulsation, mass loss and envelopes.

Figure 1 .

 1 Figure 1. Left: Number of targets observable with MATISSE in low spectral resolution (R≈35) for different types of stars, both in L band (L noft, dark blue bars) and N band (N noft, light blue bars). The transparent bars show the improvement of observability when using an external fringe tracker (L ft and N ft), like e.g. the foreseen project GRA4MAT. Right: The same plot for the high resolution of MATISSE (R≈4500 around the hydrogen Brackett α line).

Figure 2 .

 2 Figure 2. On the left, one can see modelled (Left) vs. reconstructed images (Right) of typical B[e] star disks for a 90 degrees inclination (edge-on disk). The right panel shows a 60 degrees inclination disk. The image reconstruction software used here is MiRA.

Figure 3 .

 3 Figure 3. Modelled (Left) vs. reconstructed images (Right) of a typical "pinwheel nebula" around a dusty Wolf-Rayet star. We used IRBis for this test case.

Figure 4 .

 4 Figure 4. Simulated image of the core of a R136-like cluster at the age of 2 Myr in N-band (Left) vs. reconstructed images (Right). The left panel shows 0.3" × 0.3" (FoV of MATISSE using UTs) covering 0.075 × 0.075pc 2 of the core of the cluster at the distance of LMC (50 Kpc). The image reconstruction software used here is IRBis.

Figure 5 .

 5 Figure 5. Left: the geometry of the surrounding envelop of l Car (in K band) as constrained by VINCI observations (Kervella et al. 2006). The envelope is supposed centro-symmetric and the Full-Width at Half-Maximum of the modelled gaussian is around 3 stellar radii and a flux contribution of 4%. Right: the geometry of the envelope expected in the N band is almost the same but with a larger flux contribution (figure taken from Antoine Merand PhD).

† available at http://www.ast.cam.ac.uk/ ~sverre/web/pages/nbody.htm ‡ available here: http://www.jmmc.fr

ACKNOWLEDGMENTS

This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. We used data from the Infrared Space Observatory (ISO) from ESA, and from the Infrared Astronomical Satellite (IRAS).

We are grateful to ESO, CNRS/INSU, and the Max-Planck Society for continuous support in the MATISSE project.

NN acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-15-CE31-0012-01 (project UnlockCepheids). NN acknowledges financial support from "Programme National de Physique Stellaire" (PNPS) of CNRS/INSU, France.