Sangnam Nam

Michael E Davies

Michael Elad

Rémi Gribonval

Recovery of Cosparse Signals with Greedy Analysis Pursuit in the Presence of Noise (Invited Paper)

The sparse synthesis signal model has enjoyed much success and popularity in the recent decade. Much progress ranging from clear theoretical foundations to appealing applications has been made in this field. Alongside the synthesis approach, an analysis counterpart has been used over the years. Despite the similarity, markedly different nature of the two approaches has been observed. In a recent work, the analysis model was formally formulated and the nature of the model was discussed extensively. Furthermore, a new greedy algorithm (GAP) for recovering the signals satisfying the model was proposed and its effectiveness was demonstrated.

While the understanding of the analysis model and the new algorithm has been broadened, the stability and the robustness against noise of the model and the algorithm have been mostly left out. In this work, we adapt and propose a new GAP algorithm in order to deal with the presence of noise. Empirical evidence for the algorithm is also provided.

I. INTRODUCTION

In many tasks in image and signal processing, only limited-sometimes highly insufficient-data for the signals of interest are available. Even when sufficient amount of observations is available, the data can frequently be corrupted by noise. Adding on top of these the case where the data are both incomplete and corrupted by noise, one would be tempted to simply abandon any hope of performing reasonable work on those signals.

Fortunately, the signals that we are interested in are not just 'ordinary random signals.' They exhibit highly regular behavior among them. When we mathematically model these properties shared by the signals of interest, we arrive at sometimes spectacular results even for the problems with initial bleak outlooks described above. What is an example of enabling data models?

A. The Sparse Synthesis Model

The sparse synthesis signal model is one of the most celebrated data models. Many classes of signals, e.g., images, are observed to be very low-dimensional while their actual signal dimension can be quite high. This observation is modeled as follows: Any signal x ∈ R d of interest can be expressed as a linear combination of few columns of a fixed matrix D ∈ R d×n called a dictionary. Typically, n > d so that the model is rich enough to express interesting signals.

This model has been the focus of many papers, studying its core theoretical properties by exploring practical numerical algorithms for using it in practice (e.g. [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]- [START_REF] Candès | The Dantzig selector: statistical estimation when p is much larger than n[END_REF]), evaluating theoretically these algorithms' performance guarantees, addressing ways to obtain the dictionary from a bulk of data, and beyond all these, attacking a long series of applications in signal and image processing with this model, demonstrating often state-of-the-art results.

B. The Cosparse Analysis Model

The synthesis model has a 'twin' model that takes the analysis view point. Before introducing this analysis model, let us step back and look at a typical way the sparse synthesis model is used in a signal processing application. As described at the beginning, one may have an incomplete and/or noisecontaminated observation y of a signal x. This means that there can be many-a lot of times, infinitely many-signals that could explain the observation y and hence it is difficult to recover x. However, if we know that x is sparse in D, then we can look for the coefficient vector z that is as sparse as possible while the resulting signal Dz is 'compatible' with y.

An interesting alternative has been used in practice successfully over the years; In this alternative, there is a fixed analysis operator Ω acting on the signal space, and one looks for the signal x that sparsifies the analysis representation Ωx the most while being compatible with the observation y, instead of searching for the sparsest synthesis representation z with x = Dz. Several works, e.g., [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]- [START_REF] Cai | Split bregman methods and frame based image restoration[END_REF], have empirically demonstrated the effectiveness of the analysis-based approaches and have shown either implicitly or explicitly the distinctness of two approaches (synthesis vs. analysis). On the other hand, perhaps due to the strong resemblance between the two methods, it was only recent that the underlying distinct model, which was named cosparse analysis model, was explicitly described and studied in [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF].

In the cosparse analysis model, a signal x ∈ R d is said to be cosparse with respect to an analysis operator Ω ∈ R p×d if the analysis representation Ωx contains many zero entries. Further, the number of zeros ℓ = p -Ωx 0 is called the cosparsity of x and we say that x is ℓ-cosparse. The index set where Ωx is zero will be called the cosupport of x. We want to bring to the reader's attention that we are focused on the fact that Ωx has many zeros entries rather than few non-zero entries since there is generally a non-trivial lower limit on the number of non-zero entries [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF].

C. The Greedy Analysis Pursuit Algorithm

Given an incomplete observation y of a cosparse signal x, how to recover the original signal x?

For this task, a greedy algorithm in the spirit of Matching Pursuit algorithms for the sparse synthesis model was proposed in [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF] for the cosparse signal recovery and was named Greedy Analysis Pursuit (GAP). Further details are left to Section II and [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF]. Let us now elaborate on the cosparse signal recovery problem for which the GAP algorithm was applied.

Let x ∈ R d be a cosparse signal with respect to an analysis operator Ω ∈ R p×d . Suppose that the only information available for x, apart from that it is cosparse with respect to Ω, is given by

y = Mx, (1)
where M ∈ R m×d is a measurement matrix. We assume that m < d, so the relation is underdetermined. Then, the cosparse signal recovery problem is the triplet (y, M, Ω)-possibly the quartet (y, M, Ω, ℓ), where ℓ is the coparsity of x. The goal in the cosparse signal recovery is to recover x.

The effectiveness of the GAP for the cosparse signal recovery was well-demonstrated in [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF]. Surprsingly, the GAP was observed to outperform the analysis ℓ 1 -minimization algorithm.

D. Contribution

The type of cosparse recovery problem described above rarely occurs in practice. One issue is that the clean measurements/observations y = Mx of a cosparse signal x are unlikely to be available. For this reason, it is better to suppose that y = Mx + n where n is a noise term.

Another problem of practical nature comes from the inexactness of the model. That is, most signals x of interest are actually not exactly cosparse. This leads us to work with the following approximate cosparse signal model: A signal x is approximately cosparse with respect to Ω if there is an exact cosparse signal x * close to x, i.e.,

x ≈ x * and p -Ωx * 0 ≫ 1.

Due to these practical issues, though it has delivered encouraging empirical results, the GAP algorithm has not been a viable tool to be used directly in applications. The contribution of this paper is to adapt the GAP and make it able to perform stable cosparse signal reconstruction when the observation is contaminated by noise. This is empirically demonstrated on toy MRI problems.

II. NOISELESS GAP ALGORITHM

The Greedy Analysis Pursuit (GAP) algorithm as appeared in [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF] is presented in this section.

To begin, we fix some notations. For an index set I ⊂ 1, p , Ω I is the submatrix of Ω that is obtained by retaining the rows of Ω indexed by I. A vector x 0 ∈ R d is a cosparse signal with respect to Ω and will serve as the ground truth. The cosupport of x 0 is denoted by Λ.

The GAP attempts to solve the linear inverse problem [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] where y := Mx 0 . Its aim is to identify the cosupport Λ and hence to pin down the signal x 0 . An interesting twist is that it tries to remove rows of Ω in a greedy way to arrive at Ω Λ in the end. At the first step, the GAP takes Λ0 = 1, p to be the initial estimate of the cosupport Λ and removes elements from Λ0 to obtain a new estimate Λ1 . What would be a good way to choose the elements to be discarded? If we pretend that we know x 0 , then this would be an easy task: We will compute Ω Λ0 x 0 and choose any index that corresponds to a non-zero entry. Of course, we do not have x 0 but only y. These arguments naturally suggest: First, obtain a 'good enough' estimate x0 of x 0 using y. Then, compute Ω Λ0 x0 and decide which elements to remove. For the first half of the task, we let

x0 := argmin x Ω Λ0 x 2 subject to y = Mx.
(

) 2
For the second half, note first that most likely all elements of Ω Λ0 x0 are non-zero. Hence, some care must be taken in choosing elements to remove. Assuming that x0 is a good estimate of x 0 , a safe route would be to choose the entries of Ω Λ0 x0 with the maximum magnitudes. We can see that this is a sensible approach once we view that Ω Λ0 x0 is a 'small' perturbation of Ω Λ0 x 0 . Thus, for the second part, we let

Λ1 := Λ0 \{argmax{|ω j x0 | : j ∈ Λ0 }}, (3)
where ω j is the j-th row of Ω. The version of GAP alternates the two steps above until an appropriate stopping criterion is met.

III. NEW NOISE-AWARE GAP ALGORITHM

The new version of GAP attempts to solve the linear inverse problem with noise,

y = Mx + n (4)
where y := Mx 0 + n is the noisy measurements/observations of x 0 . We assume the knowledge of the noise level ǫ; to be precise, n 2 ≤ ǫ.

Looking at how (1) is related to the estimate update step (2) of the original GAP, it is natural to consider the alternative

x0 := argmin x Ω Λ0 x 2 subject to y -Mx 2 ≤ ǫ (5)
in order to adapt the GAP for the problem (4). The same cosupport estimate update step (3) can then be used. This leads to the new noise-aware version of GAP, called GAPn. The description of the algorithm is given in Algorithm 1. a) Stopping Criteria: Various stopping criteria can be used to terminate the GAPn iterations. One criterion that proved to be useful was Ω Λi xi ∞ < η coef where η coef is some prescribed (small) value. Another possibility is xixi-1 2 < η stable for some η stable > 0. One can also decide on a prescribed cosparsity as a stopping criterion: | Λi | < ℓ stop . The first option was used in our experiment.

x′ 0 := argmin x y -Mx 2 2 + λ Ω Λ0 x 2 2
for appropriate value λ > 0. In fact, the value of λ was varied until x′ 0 satisfies the constraint y -Mx 2 ≤ ǫ. c) Super-greedy GAP Algorithms: Multiple rows of Ω can be removed at each cosupport update step. Two simple alternatives to (3) can be conceived. One way would be to discard the rows corresponding to r largest analysis coefficients at each iteration for some fixed r > 1. The second alternative is to perform

Γ i := {j ∈ Λi : |α j | ≥ t max s∈ Λi |α s |}, Λi+1 := Λi \Γ i ,
where α := Ω Λi xi , for some fixed 0 < t < 1. The Noiseless GAP using this alternative was discussed in [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF].

IV. EXPERIMENTS

We show the experimental results obtained by applying the GAPn algorithm presented in Section III.

A. Synthetic Cosparse Signal Recovery

As a first demonstration, we conducted a synthetic experiment on cosparse signal recovery. In this experiment, we fixed the signal dimension d = 200 and the size of the analysis operator p = 240. Two values of measurement size were tried: m = 160 and m = 80. The cosparsity was also fixed at ℓ = 180 and ℓ = 195 respectively for the cases m = 160 and m = 80. Once d, p, m and ℓ are fixed, the noise intensity ǫ exp were varied. At each noise intensity, the experiment was repeated 50 times and the average relative error was computed.

The details of how the cosparse signal recovery problem (y, M, Ω, ǫ) (noise intensity is added here.) was constructed are as follows: The entries of M ∈ R m×d are simply drawn independently from the normal distribution. Ω T ∈ R d×p was constructed to be a random tight frame with uniform column norm. Random ℓ elements of 1, p were selected to be a set Λ. A Gaussian random signal in the null space of Ω Λ was constructed to be x 0 . A noise term n with n 2 = ǫ exp was constructed and finally the observation y was computed by y = Mx 0 + n. The noise level ǫ := 1.01ǫ exp was passed to the algorithms. The results for the first case is shown in Figure 1. In addition to the GAPn, the analysis ℓ 1 which solves

x := argmin x Ωx 1 subject to y -Mx 2 ≤ ǫ,
was tested. In the figure, the label L1 corresponds to the result when ǫ was taken to be 0 above and the label L1n corresponds to the one where ǫ was set at the value described above. GAP denotes the noiseless version. The figure shows that while the GAP performs worse than the analysis ℓ 1 , the noise-aware GAPn outperforms the noise-aware ℓ 1 in the considered SNR range. The GAPn was also observed to work better than the ℓ 1 in the setting d = 200, p = 240, m = 160, and ℓ = 195, but the result is not shown here in the interest of space.

B. Shepp-Logan Phantom Reconstruction

As a more realistic test, we applied the GAPn for the Shepp Logan phantom reconstruction problem. In this experiment, the 256 × 256 phantom image was chosen and the measurements were made along 22 radial lines in the 2D Fourier domain. The observation was further corrupted by Gaussian noise; The noise level was ǫ exp = 0.05 Mx 0 2 and ǫ = 1.01ǫ exp , using the notations of the previous subsection.

We have tested the following algorithms: tvdantzig logbarrier and tvqc logbarrier from l1magic-1.1 package, and the GAPn. For the GAPn, the analysis operator used was the one-step finite difference operator Ω DIF that consists of all the horizontal and vertical differences of neighboring pixels.

Figure 2 demonstrates that the GAPn result captures the piece-wise constant nature of the original image better than tvdantzig or tvqc results.

C. MRI Image Reconstruction

The next experiment is related to the issue of approximate signal model. Namely, the signals we meet in practice are not exactly cosparse. For this experiment, we have chosen a (crop of) MRI image generated from the FSL MNI152 T1 0.5mm image data. 1 The setting of the problem is similar to the previous subsection, and the measurements of the image are obtained along 50 radial lines in the Fourier domain. However, no intentional noise was added.

1) Cosparsity with One-step Finite Difference Operator:

The image is not seen to be truely cosparse with respect to Ω DIF . Given this, we can envision two possible ways to reconstruct the target image. One approach is to treat the inexactness of the signal to the model as noise. Here, we start with an implicit assumption that there is a signal x * which is cosparse with respect to Ω DIF and is nearby the target signal x. Then, we recognize the problem as solving (4) with noise term given by n = M(xx * ), and try to recover x * . This requires the 'noise' intensity to be estimated.

The other strategy is to keep the constraint (1) since the observation y is thought to be clean and to apply the GAP with noise level 0 (or very small value). The algorithm stops when Ω Λi xi ∞ is deemed to be relatively small.

The resconstruction results showed that the total-variation minimization algorithms of l1magic-1.1 give better output than the GAP algorithm with analysis operator Ω DIF (result not displayed). Visually, all those results suffer from staircasing artifacts related to either 'TV-norm' or the one-step finite difference operator. The reconstruction from the GAP showed a bit more pronounced blockiness, which indicates that the GAP is enforcing the cosparse analysis model more strongly. The differences between the two approaches described above for handling the model mismatch were not significant.

2) Cosparsity with Second-order Differences: Interestingly, when the one-step difference operator was replaced by the operator Ω D2 consisting of second-order differences in the horizontal, vertical, diagonal, and anti-diagonal directions, the reconstruction result from the GAPn with ǫ = 10 -3 y 2 and η coef = 5 • 10 -4 was the best. Note that the signal x was normalized to take values in the interval [0, 1]. The improvement on SNR was from 21.0dB (Ω DIF ) to 29.6dB (Ω D2 ) which was also better than 28.4dB of the TV-based reconstruction. Perhaps more importantly, visual artifacts were significantly less pronounced. This shows the importance of using an appropriate analysis operator to fit the signals to the cosparse model. The original, and the reconstructed images using the TV-minimization and the GAPn with Ω D2 are shown in Figure 3.

V. CONCLUSION In this work, we have extended the Greedy Analysis Pursuit algorithm of [START_REF] Nam | The Cosparse Analysis Model and Algorithms[END_REF] so that it can be used for the reconstruction of cosparse signals when the noise term is present or for the reconstruction of signals that are approximately cosparse, in the context of linear inverse problems. The effectiveness of the proposed algorithm was demonstrated by experiments. The experimental results show that the GAPn algorithm outperforms the analysis ℓ 1 -minimization given the same knowledge on the noise intensity when the cosparse analysis model fits the signals of interest well. However, they also exposed the natural weakness of the algorithm when the model is not a good fit.

The weakness calls for the need to better design or learn analysis operators for signals of certain class. This need was confirmed by some improvements observed in our experimental results. The encouraging performance in this work prompts us to launch further studies: theoretical investigation of the behaviors related to the cosparse model and the algorithms; fine tuning the GAPn for the approximate cosparse model; and extensive comparisons with state-of-the-art sparse recovery algorithms (e.g., TVAL3 or ISD [START_REF] Wang | Sparse signal reconstruction via iterative support detection[END_REF]). ACKNOWLEDGMENT

Algorithm 1 :

 1 Noise-aware Greedy Analysis Pursuit (GAPn) Input: y, M, Ω, ǫ Set i = 0, Λ0 = 1, p repeat xi := argmin x Ω Λi x 2 s.t. y -Mx 2 ≤ ǫ Λi+1 := Λi \{argmax{|ω j xi | : j ∈ Λi }} i := i + 1 until stopping criteria met b) Solving the Inner Optimization (5): In the implementation of the GAPn, we solved the optimization (5) by solving an alternative problem (with conjugate gradient method):

Fig. 1 .

 1 Fig. 1. Output SNR 20 log 10 ( x 0 2 / x-x 0 2 ) as a function of input SNR 20 log 10 ( Mx 0 2 /ǫexp) for the tested algorithms for the case d = 200, p = 240, m = 160, and ℓ = 180.

Fig. 2 .

 2 Fig. 2. Shepp Logan phantom reconstruction results from 22 radial observation lines with 5% noise corruption. The original phantom (top left), tvdantzig (top right), tvqc (bottom left), and the GAP (bottom right).

Fig. 3 .

 3 Fig. 3. MRI reconstruction results from 50 radial observation lines. The original (left), TV-min. (center), and GAPn (right) reconstructed images.

Available from Debian Linux package fsl-mni152-templates.

This work was supported in part by the EU through the project SMALL (Sparse Models, Algorithms and Learning for Large-Scale data), FET-Open programme, under grant number: 225913.


Pierre-Olivier Amblard

Jean-François Coeurjolly

Frédéric Lavancier
email: frederic.lavancier@univ-nantes.fr

Anne Philippe
email: anne.philippe@univ-nantes.fr

BASIC PROPERTIES OF THE MULTIVARIATE FRACTIONAL BROWNIAN MOTION

Keywords: 1991 Mathematics Subject Classification. 26A16, 28A80, 42C40 Self similarity, Multivariate process, Long-range dependence, Superlinear process, Increment process, Limit theorem

This paper reviews and extends some recent results on the multivariate fractional Brownian motion (mfBm) and its increment process. A characterization of the mfBm through its covariance function is obtained. Similarly, the correlation and spectral analyses of the increments are investigated. On the other hand we show that (almost) all mfBm's may be reached as the limit of partial sums of (super)linear processes. Finally, an algorithm to perfectly simulate the mfBm is presented and illustrated by some simulations.

Introduction

The fractional Brownian motion is the unique Gaussian self-similar process with stationary increments. In the seminal paper of Mandelbrot and Van Ness [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF], many properties of the fBm and its increments are developed (see also [START_REF] Taqqu | Fractional Brownian motion and long-range dependence In[END_REF] for a review of the basic properties). Depending on the scaling factor (called Hurst parameter), the increment process may exhibit long-range dependence, and is commonly used in modeling physical phenomena. However in many fields of applications (e.g. neuroscience, economy, sociology, physics, etc), multivariate measurements are performed and they involve specific properties such as fractality, long-range dependence, self-similarity, etc. Examples can be found in economic time series (see [START_REF] Davidson | Alternative frequency and time domain versions of fractional Brownian motion[END_REF], [START_REF] Fleming | Wavelet-based detection of coherent structures and self-affinity in financial data[END_REF], [START_REF] Gil-Alana | A fractional multivariate long memory model for the US and the Canadian real output[END_REF]), genetic sequences [START_REF] Arianos | Cross-correlation of long-range correlated series[END_REF], multipoint velocity measurements in turbulence, functional Magnetic Resonance Imaging of several regions of the brain [START_REF] Achard | Fractal connectivity of long-memory networks[END_REF].

It seems therefore natural to extend the fBm to a multivariate framework. Recently, this question has been investigated in [START_REF] Lavancier | Covariance function of vector self-similar process[END_REF][START_REF] Lavancier | A two-sample test for comparison of long memory parameters[END_REF][START_REF] Coeurjolly | Wavelet analysis of the multivariate fractional Brownian motion[END_REF]. The aim of this paper is to summarize and to complete some of these advances on the multivariate fractional Brownian motion and its increments. A multivariate extension of the fractional Brownian motion can be stated as follows :

Definition 1. A Multivariate fractional Brownian motion (p-mfBm or mfBm) with parameter H ∈ (0, 1) p is a p-multivariate process satisfying the three following properties

• Gaussianity,

• Self-similarity with parameter H ∈ (0, 1) p ,

• Stationarity of the increments.

Here, self-similarity has to be understood as joint self-similarity. More formally, we use the following definition. Definition 2. A multivariate process (X(t) = (X 1 (t), • • • , X p (t))) t∈R is said self-similar if there exists a vector H = (H 1 , • • • , H p ) ∈ (0, 1) p such that for any λ > 0,

(1)

(X 1 (λt), • • • , X p (λt))) t∈R fidi = λ H 1 X 1 (t), • • • , λ Hp X p (t) t∈R ,
where fidi = denotes the equality of finite-dimensional distributions. The parameter H is called the self-similarity parameter.

This definition can be viewed as a particular case of operator self-similar processes by taking diagonal operators (see [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF][START_REF] Hudson | Operator-self-similar processes in a finite-dimensional space[END_REF][START_REF] Laha | Operator self-similar stochastic processes in R d[END_REF][START_REF] Maejima | Operator-self-similar stable processes[END_REF]).

Note that, as in the univariate case [START_REF] Lamperti | Semi-Stable Stochastic Processes[END_REF], the Lamperti transformation induces an isometry between the self-similar and the stationary multivariate processes. Indeed, from Definition 2, it is not difficult to check that (Y (t)) t∈R is a p-multivariate stationary process if and only if there exists H ∈ (0, 1) p such that its Lamperti transformation (t H 1 Y 1 (log(t)), . . . , t Hp Y p (log(t))) t∈R is a p-multivariate self-similar process.

The paper is organized as follows. In Section 2, we study the covariance structure of the mfBm and its increments. The cross-covariance and the cross-spectral density of the increments lead to interesting long-memory type properties. Section 3 contains the time domain as well as the spectral domain stochastic integral representations of the mfBm. Thanks to these results we obtain a characterization of the mfBm through its covariance matrix function. Section 4 is devoted to limit theorems, the mfBm is obtained as the limit of partial sums of linear processes. Finally, we discuss in Section 5 the problem of simulating sample paths of the mfBm. We propose to use the Wood and Chan's algorithm [START_REF] Wood | Simulation of stationary gaussian processes in [0, 1] d[END_REF] well adapted to generate multivariate stationary Gaussian random fields with prescribed covariance matrix function.

2. Dependence structure of the mfBm and of its increments 2.1. Covariance function of the mfBm. In this part, we present the form of the covariance matrix of the mfBm.

Firstly, as each component is a fractional brownian motion, the covariance function of the i-th component is well-known and we have

(2) EX i (s)X i (t) = σ 2 i 2 |s| 2H i + |t| 2H i -|t -s| 2H i .
with σ 2 i := var(X i (1)). The cross covariances are given in the following proposition. Proposition 3 (Lavancier et al. [START_REF] Lavancier | Covariance function of vector self-similar process[END_REF]). The cross covariances of the mfBm satisfy the following representation, for all (i, j) ∈ {1, . . . , p} 2 , i = j,

(1) If H i +H j = 1, there exists (ρ i,j , η i,j ) ∈ [-1, 1]×R with ρ i,j = ρ j,i = corr(X i (1), X j (1))
and η i,j = -η j,i such that

(3) EX i (s)X j (t) = σ i σ j 2 (ρ i,j + η i,j sign(s))|s| H i +H j + (ρ i,j -η i,j sign(t))|t| H i +H j -(ρ i,j -η i,j sign(t -s))|t -s| H i +H j .
(2) If H i +H j = 1, there exists (ρ i,j , ηi,j ) ∈ [-1, 1]×R with ρi,j = ρj,i = corr(X i (1), X j (1)) and ηi,j = -η j,i such that

(4) EX i (s)X j (t) = σ i σ j 2 {ρ i,j (|s| + |t| -|s -t|) + ηi,j (t log |t| -s log |s| -(t -s) log |t -s|)} .
Proof. Under some conditions of regularity, Lavancier et al. [START_REF] Lavancier | Covariance function of vector self-similar process[END_REF] actually prove that Proposition 3 is true for any L 2 self-similar multivariate process with stationary increments. The form of cross covariances is obtained as the unique solution of a functional equation. Formulae ( 3) and ( 4) correspond to expressions given in [START_REF] Lavancier | Covariance function of vector self-similar process[END_REF] after the following reparameterization : ρ i,j = (c i,j + c j,i )/2 and η i,j = (c i,jc j,i )/2 where c i,j and c j,i arise in [START_REF] Lavancier | Covariance function of vector self-similar process[END_REF].

Remark 1. Extending the definition of parameters ρ i,j , ρi,j , η i,j , ηi,j to the case i = j, we have ρ i,i = ρi,i = 1 and η i,i = ηi,i = 0, so that (2) coincides with (3) and (4).

Remark 2. The constraints on coefficients ρ i,j , ρi,j , η i,j , ηi,j are necessary but not sufficient conditions to ensure that the functions defined by (3) and (4) are covariance functions. This problem will be discussed in Section 3.4.

Remark 3. Note that coefficients ρ i,j , ρi,j , η i,j , ηi,j depend on the parameters (H i , H j ). Assuming the continuity of the cross covariances function with respect to the parameters (H i , H j ), the expression (4) can be deduced from (3) by taking the limit as

H i + H j tends to 1, noting that ((s + 1) H -s H -1)/(1 -H) → s log |s| -(s + 1) log |s + 1| as H → 1.
We obtain the following relations between the coefficients : as

H i + H j → 1 ρ i,j ∼ ρi,j and (1 -H i -H j )η i,j ∼ ηi,j .
This convergence result can suggest a reparameterization of coefficients η i,j in (1

-H i - H j )η i,j .
2.2. The increments process. This part aims at exploring the covariance structure of the increments of size δ of a multivariate fractional Brownian motion given by Definition 1.

Let ∆ δ X = (X(t+δ)-X(t)) t∈R denotes the increment process of the multivariate fractional Brownian motion of size δ and let ∆ δ X i be its i-th component. Let γ i,j (h, δ) = E∆ δ X i (t)∆ δ X j (t + h) denotes the cross-covariance of the increments of size δ of the components i and j. Let us introduce the function w i,j (h) given by [START_REF] Coeurjolly | Wavelet analysis of the multivariate fractional Brownian motion[END_REF] w i,j (h) = (ρ i,jη i,j sign(h))|h|

H i +H j if H i + H j = 1, ρi,j |h| + ηi,j h log |h| if H i + H j = 1.
Then from Proposition 3, we deduce that γ i,j (h, δ) is given by ( 6)

γ i,j (h, δ) = σ i σ j 2 w i,j (h -δ) -2w i,j (h) + w i,j (h + δ) .
Now, let us present the asymptotic behaviour of the cross-covariance function.

Proposition 4. As |h| → +∞, we have for any δ > 0

(7) γ i,j (h, δ) ∼ σ i σ j δ 2 |h| H i +H j -2 κ i,j (sign(h)), with (8)
κ i,j (sign(h)) = (ρ i,j -η i,j sign(h))(H i + H j )(H i + H j -1) if H i + H j = 1, ηi,j sign(h) if H i + H j = 1.
Proof. Let α = H i +H j . Let us choose h, such that |h| ≥ δ, which ensures that sign(h-δ) = sign(h) = sign(h + δ). When α = 1, this allows us to write

γ i,j (h, δ) = σ i σ j 2 |h| α (ρ i,j -η i,j sign(h)) B(h), with B(h) = 1 -δ h α -2 + 1 + δ h α ∼ α(α -1)δ 2 h -2 , as |h| → +∞. When α = 1 and |h| ≥ δ, γ i,j (h, δ) reduces to γ i,j (h, δ) = σ i σ j 2 ηi,j B(h) with B(h) = (h -δ) log 1 - δ h + (h + δ) log 1 + δ h .
Using the expansion of log(1 ± x) as x → 0 leads to B(h) ∼ δ 2 |h| -1 as |h| → +∞, which implies the result.

Proposition 4 and (6) lead to the following important remarks on the dependence structure. For i = j and H i + H j = 1 :

• If the two fractional Gaussian noises are short-range dependent (i.e. H i < 1/2 and H j < 1/2) then they are either short-range interdependent if ρ i,j = 0 or η i,j = 0, or independent if ρ i,j = η i,j = 0. • If the two fractional Gaussian noises are long-range dependent (i.e. H i > 1/2 and H j > 1/2) then they are either long-range interdependent if ρ i,j = 0 or η i,j = 0, or independent if ρ i,j = η i,j = 0. This confirms the dichotomy principle observed in [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF]. • In the other cases, the two fractional Gaussian noises can be short-range interdependent if ρ i,j = 0 or η i,j = 0 and H i + H j < 1, long-range interdependent if ρ i,j = 0 or η i,j = 0 and H i + H j > 1 or independent if ρ i,j = η i,j = 0.

Moreover, note that when H i + H j = 1, whatever the nature of the two fractional Gaussian noises (i.e. short-range or long-range dependent, or even independent), they are either long-range interdependent if ηi,j = 0 or independent if ηi,j = 0.

The following result characterizes the spectral nature of the increments of a mfBm.

Proposition 5 (Coeurjolly et al. [START_REF] Coeurjolly | Wavelet analysis of the multivariate fractional Brownian motion[END_REF]). Let S i,j (•, δ) be the (cross)-spectral density of the increments of size δ of the components i and j, i.e. the Fourier transform of γ i,j (•, δ) S i,j (ω, δ) = 1 2π R e -ihω γ i,j (h, δ) dh =: F T (γ i,j (•, δ)).

(i) For all i, j and for all H i , H j , we have

(9) S i,j (ω, δ) = σ i σ j π Γ(H i + H j + 1) 1 -cos(ωδ) |ω| H i +H j +1 × τ i,j (sign(ω)), where (10)
τ i,j (sign(ω)) = ρ i,j sin π 2 (H i + H j ) -iη i,j sign(ω) cos π 2 (H i + H j ) if H i + H j = 1, ρi,j -i π 2 ηi,j sign(ω) if H i + H j = 1.
(ii) For any fixed δ, when H i + H j = 1 then we have, as ω → 0,

S i,j (ω, δ) ∼ σ i σ j 2π Γ(H i + H j + 1)δ 2 ρ 2 i,j sin π 2 (H i + H j ) 2 + η 2 i,j cos π 2 (H i + H j ) 2 1/2 |ω| H i +H j -1 . (11)
(iii) Moreover, when H i + H j = 1, the coherence function between the two components i and j satisfies, for all ω

C i,j (ω, δ) := |S i,j (ω, δ)| 2 S i,i (ω, δ)S j,j (ω, δ) = Γ(H i + H j + 1) 2 Γ(2H i + 1)Γ(2H j + 1)
ρ 2 i,j sin π 2 (H i + H j ) 2 + η 2 i,j cos π 2 (H i + H j ) 2 sin(πH i ) sin(πH j ) . ( 12
) (iv) When H i +H j = 1, (11) and (12) hold, replacing ρ 2 i,j sin π 2 (H i + H j ) 2 +η 2 i,j cos π 2 (H i + H j ) 2 by ρ2 i,j + π 2 4 η2 i,j . Proof.
The proof is essentially based on the fact that in the generalized function sense, for α > -1,

F T (|h| α ) = - 1 π Γ(α + 1) sin π 2 α |ω| -α-1 , F T (h α + ) = 1 2π Γ(α + 1)e -isign(ω) π 2 (α+1) |ω| -α-1 , F T (h α -) = 1 2π Γ(α + 1)e isign(ω) π 2 (α+1) |ω| -α-1 , F T (h log |h|) = i sign(ω) 2ω 2 .
See [START_REF] Coeurjolly | Wavelet analysis of the multivariate fractional Brownian motion[END_REF] for more details.

Remark 4. From this proposition, we retrieve the same properties of dependence and interdependence of X i and X j as stated after Proposition 4.

Time reversibility.

A stochastic process is said to be time reversible if X(t) = X(-t) for all t. As shows in [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF], this is equivalent for zero-mean multivariate Gaussian stationary processes to EX i (t)X j (s) = EX i (s)X j (t) for s, t ∈ R or that the cross covariance of the increments satisfies γ i,j (h, δ) = γ i,j (-h, δ) for h ∈ R. The following proposition characterizes this property. Proposition 6. A mfBm is time reversible if and only if η i,j = 0 (or ηi,j = 0) for all i, j = 1, . . . , p.

Proof. If η i,j = 0 (or ηi,j = 0), γ i,j (h, δ) is proportional to the covariance of a fractional Gaussian noise with Hurst parameter (H i + H j )/2 and is therefore symmetric. Let us prove the converse. Let α = H i + H j , then

γ i,j (h, δ) -γ i,j (-h, δ) = σ i σ j × -η i,j (sign(h -δ)|h -δ| α + 2sign(h)|h| α -sign(h + δ)|h + δ| α ) if α = 1, ηi,j ((h -δ) log |h -δ| -2h log |h| + (h + δ) log |h + δ|) if α = 1.
Assuming γ i,j (h, δ)γ i,j (-h, δ) equals zero for all h leads to η i,j = 0 (or ηi,j = 0).

Remark 5. This result can also be viewed from a spectral point view. The time reversibility of a mfBm is equivalent to the fact that the spectral density matrix is real. Using (9), this implies η i,j = 0 (or ηi,j = 0).

Integral representation

3.1. Spectral representation. The following proposition contains the spectral representation of mfBm. This representation will be especially useful to obtain a condition easy to verify which ensures that the functions defined by ( 3) and ( 4) are covariance functions.

Theorem 7 (Didier and Pipiras, [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF]). Let (X(t)) t∈R be a mfBm with parameter (H 1 , • • • , H p ) ∈ (0, 1) p . Then there exists a p × p complex matrix A such that each component admits the following representation

(13) X i (t) = p j=1 e itx -1 ix (A ij x -H i +1/2 + + Āij x -H i +1/2 - ) Bj ( dx),
where for all j = 1, . . . , p, Bj is a Gaussian complex measure such that Bj = Bj,1 + i Bj,2 with Bj,1 (x) = Bj,1 (-x), Bj,2 (x) = -Bj,2 (x), Bj,1 and Bj,2 are independent and

E( Bj,i ( dx) Bj,i ( dx) ′ ) = dx, i = 1, 2.
Conversely, any p-multivariate process satisfying (13) is a mfBm process.

Proof. This representation is deduced from the general spectral representation of operator fractional Brownian motions obtained in [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF]. By denoting

-H + 1/2 := diag(-H 1 + 1/2, • • • , -H p + 1/2) we have indeed (14) X(t) = e itx -1 ix (x -H+1/2 + A + x -H+1/2 - Ā) B( dx).
Any mfBm having representation (13) has a covariance function as in Proposition 3. The coefficients ρ i,j , η i,j , ρi,j and ηi,j involved in (3) and ( 4) satisfy [START_REF] Gil-Alana | A fractional multivariate long memory model for the US and the Canadian real output[END_REF] (AA

* ) i,j = σ i σ j 2π Γ(H i + H j + 1)τ i,j (1)
,

where τ i,j is given in [START_REF] Davidson | The functional central limit theorem and weak convergence to stochastic integrals[END_REF] and where A * is the transpose matrix of Ā. This relation is obtained by identification of the spectral matrix of the increments deduced on the one hand from ( 13) and provided on the other hand in Proposition 5. Given [START_REF] Dietrich | Fast and exact simulation of stationary gaussian processes through circulant embedding of the covariance matrix[END_REF], relation [START_REF] Gil-Alana | A fractional multivariate long memory model for the US and the Canadian real output[END_REF] provides easily the coefficients ρ i,j , η i,j , ρi,j and ηi,j which define the covariance function. The converse is more difficult to obtain. Given a covariance function as in Proposition 3, obtaining the explicit representation (13) requires finding a matrix A such that (15) holds. This choice is possible if and only if the matrix on the right hand side of ( 15) is positive semidefinite. Then a matrix A (which is not unique) may be deduced by the Cholesky decomposition. When p = 2, an explicit solution is the matrix with entries, for i, j = 1, 2,

A i,j = λ i,j ρ i,j sin π 2 (H i + H j ) + η i,j 1 -C i,j C i,j cos π 2 (H i + H j ) +i ρ i,j 1 -C i,j C i,j sin π 2 (H i + H j ) -η i,j cos π 2 (H i + H j ) ,
where

λ i,j = σ i 2 √ π Γ(H i + H j + 1) Γ(2H j + 1) sin(H j π)
and C i,j is given in [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF], provided

H 1 + H 2 = 1.
When H 1 + H 2 = 1, the same solution holds, replacing ρ i,j by ρi,j and η i,j cos π 2 (H 1 + H 2 ) by -π 2 ηi,j .

Moving average representation.

In the next proposition, we give an alternative characterization of the mfBm from an integral representation in the time domain (or moving average representation).

Theorem 8 (Didier and Pipiras, [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF]). Let (X(t)) t∈R be a mfBm with parameter (H 1 , • • • , H p ) ∈ (0, 1) p . Assume that for all i ∈ {1, ..., p}, H i = 1/2. Then there exist M + , M -two p × p real matrices such that each component admits the following representation (16)

X i (t) = p j=1 R M + i,j (t -x) H i -.5 + -(-x) H i -.5 + + M - i,j (t -x) H i -.5 - -(-x) H i -.5 - W j (dx), with W (dx) = (W 1 (dx), • • • , W p (dx)
) is a Gaussian white noise with zero mean, independent components and covariance EW i (dx)W j (dx) = δ i,j dx. Conversely, any p-multivariate process satisfying (16) is a mfBm process.

Proof. This representation is deduced from the general representation obtained in [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF]. Remark 6. When H i = 1/2 for each i ∈ {1, ..., p}, it is shown in [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF] that each component of the mfBm admits the following representation :

X i (t) = p j=1 R M + i,j (sign(t -x) -sign(x)) + M - i,j (log |t -x| -log |x|) W j (dx).
Our conjecture is that this representation remains valid when H i = 1/2 whatever the values of other parameters H j , j = i.

Starting from the moving average representation ( 16), using results in [START_REF] Stoev | How rich is the class of multifractional Brownian motions? Stochastic Processes and their Applications[END_REF], we can specify the coefficients ρ i,j , η i,j , ρi,j and ηi,j involved in the covariances (3) and (4) (see [START_REF] Lavancier | Covariance function of vector self-similar process[END_REF]). More precisely, let us denote

M + (M + ) ′ = α ++ i,j , M -(M -) ′ = α -- i,j , M + (M -) ′ = α +- i,j
where M ′ is the transpose matrix of M. The variance of each component is equal to

σ 2 i = B(H i + .5, H i + .5) sin(H i π) α ++ i,i + α -- i,i -2 sin(H i π)α +- i,i
, where B(•, •) denotes the Beta function.

Moreover, if H i + H j = 1 then

σ i σ j ρ i,j = B(H i + .5, H j + .5) sin((H i + H j )π) × (α ++ i,j + α -- i,j )(cos(H i π) + cos(H j π)) -(α +- i,j + α -+ i,j ) sin((H i + H j )π) , σ i σ j η i,j = B(H i + .5, H j + .5) sin((H i + H j )π) × (α ++ i,j -α -- i,j )(cos(H i π) -cos(H j π)) -(α +- i,j -α -+ i,j ) sin((H i + H j )π) . If H i + H j = 1 then σ i σ j ρi,j = B(H i + .5, H j + .5) sin(H i π) + sin(H j π) 2 (α ++ i,j + α -- i,j ) -α +- i,j -α -+ i,j , σ i σ j ηi,j = (H j -H i )(α ++ i,j -α -- i,j
). Conversely, given a covariance function as in Proposition 3, if H i = 1/2 for all i, one may find matrices M + and M -such that (16) holds, provided the matrix on the right hand side of ( 15) is positive semidefinite. Indeed, in this case, a matrix A which solves (15) may be found by the Cholesky decomposition, then M + and M -are deduced from relation (3.20) in [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF]:

M ± = π 2 D -1 1 A 1 ± D -1 2 A 2 ,
where A = A 1 + iA 2 and

D 1 = diag sin(πH 1 )Γ(H 1 + 1 2 ), . . . , sin(πH p )Γ(H p + 1 2 ) , D 2 = diag cos(πH 1 )Γ(H 1 + 1 2 ), . . . , cos(πH p )Γ(H p + 1 2
) .

3.3.

Two particular examples. Let us focus on two particular examples which are quite natural: the causal mfBm (M -= 0) and the well-balanced mfBm (M -= M + ). In the causal case, the integral representation is a direct generalization of the integral representation of Mandelbrot and Van Ness [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF] to the multivariate case. The well-balanced case is studied by Stoev and Taqqu in one dimension [START_REF] Stoev | How rich is the class of multifractional Brownian motions? Stochastic Processes and their Applications[END_REF]. With the notation of the two previous sections, we note that the causal case (resp. well-balanced case) leads to A 1 = tan(πH)A 2 (resp. A 2 = 0), where tan(πH) := diag (tan(πH 1 ), . . . , tan(πH p )). In these two cases, the covariance only depends on one parameter, for instance ρ i,j (or ρi,j ). Indeed we easily deduce η i,j (or ηi,j ) as follows :

• in the causal case i.e. M -= 0 or equivalently A 1 = tan(πH)A 2 :

η i,j = -ρ i,j tan( π 2 (H i + H j )) tan( π 2 (H i -H j )) if H i + H j = 1, ηi,j = ρi,j 2 π tan(πH i ) if H i + H j = 1.
• in the well-balanced case i.e. M -= M + or equivalently A 2 = 0 :

η i,j = 0 if H i + H j = 1, ηi,j = 0 if H i + H j = 1.
Remark 7. From Proposition 6, the well-balanced mfBm is therefore time reversible.

3.4.

Existence of the covariance of the mfBm. In this paragraph, we highlight some of the previous results in order to exhibit the sets of the possible parameters (ρ i,j , η i,j ) or (ρ i,j , ηi,j ) ensuring the existence of the covariance of the mfBm. For i, j = 1, . . . , p, let us give

(H i , H j ) ∈ (0, 1) 2 , (σ i , σ j ) ∈ R + × R + and (ρ i,j η i,j ) ∈ [-1, 1] × R with ρ j,i = ρ i,j and η j,i = -η i,j if H i + H j = 1, or (ρ i,j , ηi,j ) ∈ [-1, 1] × R with ρj,i = ρi,j and ηj,i = -η i,j if H i + H j = 1.
For this set of parameters, let us define the matrix Σ(s, t) = (Σ i,j (s, t)) as follows : Σ i,i (s, t) is given by (2) and Σ i,j (s, t) is given by (3) when H i + H j = 1 and (4) when

H i + H j = 1.
Proposition 9. The matrix Σ(s, t) is a covariance matrix function if and only if the Hermitian matrix Q = (Γ(H i + H j + 1)τ i,j (1)) with τ i,j defined in [START_REF] Davidson | The functional central limit theorem and weak convergence to stochastic integrals[END_REF], is positive semidefinite. When p = 2, this condition reduces to C 1,2 ≤ 1 where C 1,2 is the coherence defined by [START_REF] Didier | Integral representations of operator fractional Brownian motion[END_REF].

Proof. First, note that since ρ j,i = ρ i,j and η j,i = -η i,j , Q is a Hermitian matrix. Now, if Q is positive semidefinite, then so is the matrix (2π) -1 (σ i σ j Q i,j ). Therefore there exists a matrix A satisfying [START_REF] Gil-Alana | A fractional multivariate long memory model for the US and the Canadian real output[END_REF]. From Theorem 7, there exists a mfBm having Σ(s, t) as covariance matrix function. Conversely, if Σ(s, t) is a covariance matrix function of a mfBm then the representation (13) holds and by [START_REF] Gil-Alana | A fractional multivariate long memory model for the US and the Canadian real output[END_REF], the matrix Q is positive semidefinite.

When p = 2, the result comes from the fact that Q is positive semidefinite if and only if det(Q) ≥ 0 or equivalently C 1,2 ≤ 1.

When p = 2, for fixed values of (H 1 , H 2 ) the condition C 1,2 ≤ 1 means that the set of possible parameters (ρ 1,2 , η 1,2 ) is the interior of an ellipse. These sets are represented in Figure 1 according to different values of H 1 and H 2 . Note that, in order to compare the cases

H 1 +H 2 = 1 and H 1 +H 2 = 1, we have reparameterized η 1,2 by η ′ 1,2 := η 1,2 (1-H 1 -H 2 )
. In such a way, the second ellipse becomes the limit of the first one as H 1 + H 2 → 1 (see also Remark 3).

Let us underline that the maximum possible correlation between two fBm's is obtained when η 1,2 = 0, i.e. when the 2-mfBm is time reversible according to Proposition 6. Remark 8. When H 1 = . . . = H p = H = 1/2, the matrix Q rewrites Q i,j = Γ(2H + 1)(sin(πH)ρ i,jiη i,j cos(πH)) and • if the mfBm is time reversible, i.e. η i,j = 0 (for i, j = 1, . . . , p), then Q is a correlation matrix and is therefore positive-semidefinite for any |ρ i,j | ≤ 1, • when p = 2, the set of possible values for (ρ 1,2 , η 1,2 ) associated to H and 1 -H are the same.

In the particular case of the causal or the well-balanced mfBm, the matrix Σ(s, t) can be expressed through the sole parameter ρ i,j . Let λ(H 1 , H 2 ) the function which equals to cos( π 2 (H 1 -H 2 )) 2 in the causal case and which equals 1 in the well-balanced case. The maximal possible correlation when p = 2 is given by

ρ 2 1,2 = Γ(2H 1 + 1)Γ(2H 2 + 1) Γ(H 1 + H 2 + 1) 2 sin(πH 1 ) sin(πH 2 ) sin( π 2 (H 1 + H 2 )) 2 × λ(H 1 , H 2 ).
(ρ 1,2 , η ′ 1,2 ) with η ′ 1,2 := η 1,2 (1 -H 1 -H 2 ) when H 1 + H 2 = 1 and (ρ 1,2 , η1,2 ) when H 1 + H 2 = 1,
ensuring that Σ(s, t) is a covariance matrix function in the particular case p = 2.

The mfBm as a limiting process.

A natural way to construct self-similar processes is through limits of stochastic processes. In dimension one, the result is due to Lamperti [START_REF] Lamperti | Semi-Stable Stochastic Processes[END_REF]. In [START_REF] Hudson | Operator-self-similar processes in a finite-dimensional space[END_REF], an extension to operator selfsimilar processes is given. A similar result for the mfBm is deduced and stated below. In the following, a p-multivariate process (X(t)) t∈R is said proper if, for each t, the law of X(t) is not contained in a proper subspace of R p . Theorem 10. Let (X(t)) t∈R be a p-multivariate proper process, continuous in probability. If there exist a p-multivariate process (Y (t)) t∈R and p real functions a 1 , ...., a p such that

(17) (a 1 (n)Y 1 (nt), . . . , a p (n)Y p (nt)) n→∞ ---→ f idi X(t),
then the multivariate process (X(t)) is self-similar. Conversely, any multivariate selfsimilar process can be obtained as a such limit.

Proof. The proof is similar to Theorem 5 in [START_REF] Hudson | Operator-self-similar processes in a finite-dimensional space[END_REF]. Fix k ∈ N and r > 0. For each T ∈ R k we denote X(T ) := (X(T 1 ), . . . , X(T k )). Let D r,k be the set of all invertible diagonal matrices α such that, for all T ∈ R k , X(rT ) = αX(T ).

Let us first show that D r,k is not empty. According to [START_REF] Laha | Operator self-similar stochastic processes in R d[END_REF], we have

diag(a 1 (n), . . . , a p (n))Y (nrT ) n→∞ ---→ d X(rT ),
and diag(a 1 (rn), . . . , a p (rn))Y (nrT

) n→∞ ---→ d X(T ).
Since (X(t)) is proper, diag(a 1 (n), . . . , a p (n)) and diag(a 1 (rn), . . . , a p (rn)) are invertible for n large enough. Then, Theorem 2.3 in [START_REF] Weissman | On convergence of types and processes in Euclidean space[END_REF] ensures that α n defined by α n = diag(a 1 (n), . . . , a p (n))diag(a 1 (nr), . . . , a p (nr)) -1 has a limit in D r,k . Moreover if α is a limit of α n then X(rT ) = αX(T ) and thus D r,k = ∅.

It is then straightforward to adapt Lemma 7.2-7.5 in [START_REF] Hudson | Operator-self-similar processes in a finite-dimensional space[END_REF] for the subgroup D r,k , which yields that for each r, ∩ k D r,k is not empty. Therefore, for any fixed r > 0, there exists α ∈ ∩ k D r,k such that (X(rt)) and (αX(T )) have the same finite dimensional distributions. Theorem 1 in [START_REF] Hudson | Operator-self-similar processes in a finite-dimensional space[END_REF] ensures that there exists (H 1 , . . . , H p ) ∈ (0, 1) p such that α = diag(r H 1 , . . . , r Hp ). The converse is trivial.

As an illustration of Theorem 10, the mfBm can be obtained as the weak limit of partial sums of sum of linear processes (also called superlinear processes, see [START_REF] Zhao | On martingale approximations[END_REF]). Some examples may be found in [START_REF] Chung | Sample means, sample autocovariances, and linear regression of stationary multivariate long memory processes[END_REF] and [START_REF] Lavancier | A two-sample test for comparison of long memory parameters[END_REF]. In Proposition 11 below, we give a general convergence result which allows to reach almost any mfBm from such partial sums. The unique restriction concerns the particular case when at least one of the Hurst parameters is equal to 1/2.

Let (ǫ j (k)) k∈Z , j = 1, . . . , p be p independent i.i.d. sequences with zero mean and unit variance. Let us consider the superlinear processes ( 18)

Z i (t) = p j=1 k∈Z ψ i,j (t -k)ǫ j (k), i = 1, . . . , p,
where ψ i,j (k) are real coefficients with k∈Z ψ 2 i,j (k) < ∞. Moreover, we assume that ψ i,j (k) = ψ + i,j (k) + ψ - i,j (k) where ψ + i,j (k) satisfies one of the following conditions:

(i) ψ + i,j (k) = α + i,j + o(1) k d + i,j -1 + as |k| → ∞, with 0 < d + i,j < 1 2 and α + i,j = 0, (ii) ψ + i,j (k) = α + i,j + o(1) k d + i,j -1 +
as |k| → ∞, with -1 2 < d + i,j < 0, k∈Z ψ + i,j (k) = 0 and α + i,j = 0, (iii) k∈Z ψ + i,j (k) < ∞ and let α + i,j := k∈Z ψ + i,j (k) = 0, d + i,j := 0. Similarly, ψ - i,j (k) is assumed to satisfy (i), (ii) or (iii) where k + , d + i,j and α + i,j are replaced by k -, d - i,j and α - i,j .

Proposition 11. Let d i = max(d + i1 , d - i1 , • • • , d + ip , d - ip )
, for i = 1, . . . , p. Consider the vector of partial sums, for τ ∈ R,

S n (τ ) =   n -d 1 -(1/2) [nτ ] t=1 Z 1 (t), • • • , n -dp-(1/2) [nτ ] t=1 Z p (t)   .
Then the finite dimensional distributions of (S n (τ )) τ ∈R converge in law towards a p-mfBm (X(τ )) τ ∈R .

• When d i = 0, (X i (τ )) τ ∈R is defined through the integral representation [START_REF] Hudson | Operator-self-similar processes in a finite-dimensional space[END_REF] where

M + i,j = α + i,j d -1 i 1 d + i,j =d i and M - i,j = α - i,j d -1 i 1 d - i,j =d i . • When d i = 0, X i (τ ) = p j=1 (α + i,j 1 d + i,j =0 + α - i,j 1 d - i,j =0 )W j (τ )
, where W j is a standard Brownian motion.

Moreover, if for all j = 1, . . . , p, E(ǫ j (0) 2α ) < ∞ with α > 1 ∨ (1 + 2d max ) -1 where d max = max i {d i }, then S n (.) converges towards the p-mfBm X(.) in the Skorohod space D([0, 1]). Sketch of proof. We focus on the convergence in law of S n (τ ) to X(τ ), for a fixed τ in R, the finite dimensional convergence is deduced in the same way. We set for simplicity τ = 1.

According to the Cramér-Wold device, for any vector (λ 1 , . . . , λ p ) ∈ R p , we must show that λ ′ S n (1) converges in law to λ ′ X(1). We may rewrite λ ′ S n (1) as a sum of discrete stochastic integrals (see [START_REF] Surgailis | Non CLTs: U-statistics, multinomial formula and approximations of multiple Itô-Wiener integrals[END_REF] and [START_REF] Bružaitė | Asymptotic independence of distant partial sums of linear process[END_REF]) :

λ ′ S n (1) = p i=1 λ i n -d i -(1/2) n t=1 Z i (t) = p i=1 λ i p j=1 k∈Z n -d i -(1/2) n t=1 (ψ + i,j (t -k) + ψ - i,j (t -k))ǫ j (k) = p i=1 λ i p j=1 R (f + i,j,n (x) + f - i,j,n (x))W j,n (dx), (19)
where the stochastic measures W j,n , j = 1, . . . , p are defined on finite intervals C by

W j,n (C) = n -1/2 k/n∈C ǫ j (k),
and where f + i,j,n , f - i,j,n are piecewise constant functions defined as follows: denoting ⌈x⌉ the smallest integer not less than x, we have for all x ∈ R,

f + i,j,n (x) = n -d i n t=1 ψ + i,j (t -⌈nx⌉), respectively f - i,j,n (x) = n -d i n t=1 ψ - i,j (t -⌈xn⌉).
The following lemma states the convergence of a linear combination of discrete stochastic integrals as in [START_REF] Lavancier | A two-sample test for comparison of long memory parameters[END_REF]. A function is said n-simple if it takes a finite number a nonzero constant values on intervals (k/n, (k + 1)/n], k ∈ Z.

Lemma 12. Let (f 1,n , • • • , f p,n ) n∈N be a sequence of p n-simple functions in L 2 (R). If for any j = 1, . . . , p, there exists f j ∈ L 2 (R) such that R |f j,n (x) -f j (x)| 2 dx → 0, then p j=1 R f j,n (x)W j,n (dx) converges in law to p j=1 R f j (x)W j (dx)
, where the W j 's are independent standard Gaussian random measures.

When p = 1, this lemma is proved in [START_REF] Surgailis | Zones of attraction of self-similar multiple integrals[END_REF]. The case p = 2 is considered in [START_REF] Bružaitė | Asymptotic independence of distant partial sums of linear process[END_REF] and the extension to p ≥ 3 is straightforward.

From Lemma 12 and (19), it remains to show that

lim n→∞ R f ± i,j,n (x) - α ± i,j d i (1 -x) d i ± -(-x) d i ± 1 d ± i,j =d i 2 dx = 0,
where we agree that d

-1 i ((1-x) d i ± -(-x) d i ± ) = 1 [0,1] (x) when d i = 0.
Below, we only consider the pointwise convergence of f ± i,j,n (x), for x ∈ R, when d ± i,j = d i . The convergence in L 2 is then deduced from the dominated convergence theorem (see [START_REF] Surgailis | Zones of attraction of self-similar multiple integrals[END_REF], [START_REF] Surgailis | Non CLTs: U-statistics, multinomial formula and approximations of multiple Itô-Wiener integrals[END_REF], [START_REF] Bružaitė | Asymptotic independence of distant partial sums of linear process[END_REF] for details). It also follows easily that, when

d ± i,j < d i , R |f ± i,j,n (x)| 2 dx → 0. Under assumption (i), note that since d i > 0, (1 -x) d i ± -(-x) d i ± = d i 1 0 (t -x) d i -1 ± dt. We have, for any x ∈ R, f ± i,j,n (x) = n -d i n t=1 ψ ± i,j (t -⌈nx⌉) = n -d i n 0 ψ ± i,j (⌈t⌉ -⌈nx⌉)dt = n -d i n 0 α ± i,j + o(1) (⌈t⌉ -⌈nx⌉) d i -1 ± dt = 1 0 α ± i,j + o(1) ⌈nt⌉ -⌈nx⌉ n d i -1 ± dt -→ α ± i,j 1 0 (t -x) d i -1 ± dt.
Under assumption (ii),

d i < 0. When x ≤ 0, (1 -x) d i + -(-x) d i + = d i 1 0 (t -x) d i -1
dt and the convergence of f + i,j,n (x) can be proved as above. When x ≥ 1, (1x)

d i + -(-x) d i + = 0 = f + i,j,n (x). When 0 ≤ x ≤ 1, (1 -x) d i + -(-x) d i + = -d i +∞ 1
(tx) d i -1 dt and, since k∈Z ψ + i,j (k) = 0, we have

f + i,j,n (x) = n -d i n t=⌈nx⌉ ψ + i,j (t -⌈nx⌉) = n -d i n-⌈nx⌉ t=0 ψ + i,j (t) = -n -d i t>n-⌈nx⌉ ψ + i,j (t).
Therefore,

f + i,j,n (x) = -n -d i +∞ n-⌈nx⌉ α + i,j + o(1) (⌈t⌉) d i -1 dt = - +∞ 1-⌈nx⌉ n α + i,j + o(1) ⌈nt⌉ n d i -1 dt -→ -α + i,j +∞ 1-x t d i -1 dt = -α + i,j +∞ 1 (t -x) d i -1 dt This proves f + i,j,n (x) → d -1 i α + i,j ((1 -x) d i + -(-x) d i + )
, for any x ∈ R, under assumption (ii). The same scheme may be used to prove that f

- i,j,n (x) → d -1 i α - i,j ((1 -x) d i --(-x) d i -)
under assumption (ii), noting that

(1 -x) d i --(-x) d i -=      0 when x ≤ 0, -d i 0 -∞ (t -x) d i -1 dt when 0 ≤ x ≤ 1, d i 1 0 (t -x) d i -1 - dt when x > 1.
Under assumption (iii),

f ± i,j,n (x) = n t=1 ψ ± i,j (t -⌈nx⌉) = n-⌈nx⌉ t=1-⌈nx⌉ ψ ± i,j (t). Since t∈Z ψ ± i,j (t) < ∞, f ± i,j,n (x) → 0 for all x / ∈ [0, 1]. When x ∈ [0, 1], we have f ± i,j,n (x) → α ± i,j .
Therefore, the first claim of the theorem is proved, i.e. the convergence in law of the finite dimensional distribution of (S n (τ )) τ ∈R to (X(τ )) τ ∈R . To extend this convergence to a functional convergence in D([0, 1]), it remains to show tightness of the sequence (S n (τ )) τ ∈[0,1] . This follows exactly from the same arguments as in the proof of Theorem 1.2 in [START_REF] Bružaitė | Asymptotic independence of distant partial sums of linear process[END_REF].

5. Synthesis of the mfBm 5.1. Introduction. The exact simulation of the fractional Brownian motion has been a question of great interest in the nineties. This may be done by generating a sample path of a fractional Gaussian noise. An important step towards efficient simulation was obtained after the work of Wood and Chan [START_REF] Wood | Simulation of stationary gaussian processes in [0, 1] d[END_REF] about the simulation of arbitrary stationary Gaussian sequences with prescribed covariance function. The technique relies upon the embedding of the covariance matrix into a circulant matrix, a square root of which is easily calculated using the discrete Fourier transform. This leads to a very efficient algorithm, both in terms of computation time and storage needs. Wood and Chan method is an exact simulation method provided that the circulant matrix is semidefinite positive, a property that is not always satisfied. However, for the fractional Gaussian noise, it can be proved that the circulant matrix is definite positive for all H ∈ (0, 1), see [START_REF] Craigmile | Simulating a class of stationary gaussian processes using the davies-harte algortihm, with application to long memory processes[END_REF][START_REF] Dietrich | Fast and exact simulation of stationary gaussian processes through circulant embedding of the covariance matrix[END_REF].

In [START_REF] Chan | Simulation of stationary gaussian vector fields[END_REF], Wood and Chan extended their method and provided a more general algorithm adapted to multivariate stationary Gaussian processes. The main characteristic of this method is that if a certain condition for a familiy of Hermitian matrices holds then the algorithm is exact in principle, i.e. the simulated data have the true covariance. We present hereafter the main ideas, briefly describe the algorithm and propose some examples. Remark 9. Other approaches could have been undertaken (see [START_REF] Bardet | Generators of long-range dependent processes : A survey[END_REF] for a review in the case p = 1). Approximate simulations can be done by discretizing the moving-average or spectral stochastic integrals [START_REF] Dietrich | Fast and exact simulation of stationary gaussian processes through circulant embedding of the covariance matrix[END_REF] or [START_REF] Hudson | Operator-self-similar processes in a finite-dimensional space[END_REF]. [START_REF] Chambers | The simulation of random vector time series with given spectrum[END_REF] also proposed an approximate method based on the spectral density matrix of the increments for synthesizing multivariate Gaussian time series. Thanks to Proposition 5, this could also be envisaged for the mfBm.

Method and algorithm.

For two arbitrary matrices A = (A j,k ) and B, we use A⊗B to denote the Kronecker product of A and B that is the block matrix (A j,k B).

Let ∆X := ∆ 1 X denotes the increments of size 1 (δ = 1) of a mfBm. We have ∆X = (∆X(t)) t∈R = ((∆X 1 (t), . . . , ∆X p (t)) ′ ) t∈R . The aim is to simulate a realization of a multivariate fractional Gaussian noise discretized at times j = 1, . . . , n, that is a realization of (∆X(1), . . . , ∆X(n)). Then a realization of the discretized mfBm will be easily obtained.

We denote by ∆X (n) the merged vector ∆X (n) = (∆X(1) ′ , . . . , ∆X(n) ′ ) ′ and by G its covariance matrix. G is the np×np Toeplitz block matrix G = (G(|i-j|)) i,j=1,...,n where for h = 0, . . . , n-1, G(h) is the p×p matrix given by G(h) := (γ j,k (h)) j,k=1,...,p . The simulation problem can be viewed as the generation of a random vector following a N np (0, G). This may be done by computing G 1/2 but the complexity of such a procedure is O(p 3 n 2 ) for block Toeplitz matrices. To overcome this numerical cost, the idea is to embed G into the block circulant matrix C = circ{C(j), j = 0, . . . , m -1}, where m is a power of 2 greater than 2(n -1) and where each C(j) is the p × p matrix defined by ( 20)

C(j) =    G(j) if 0 ≤ j < m/2 1 2 G(j) + G(j) ′ if j = m/2 G(j -m) if m/2 < j ≤ m -1.
Such a definition ensures that C is a symmetric matrix with nested block circulant structure and that G = {C(j), j = 0, . . . , n -1} is a submatrix of C. Therefore, the simulation of a N np (0, G) may be achieved by taking the n "first" components of a vector N mp (0, C), which is done by computing C 1/2 . The last problem is more simple since one may exploit the circulant characteristic of C: there exist m Hermitian matrices B(j) of size p × p such that the following decomposition holds

(21) C = (J ⊗ I p ) diag(B(j), j = 0, . . . , m -1) (J * ⊗ I p ),
where Q is the m × m unitary matrix defined for j, k = 0, m -1 by J j,k = e -2iπjk/m . The computation of C 1/2 is much less expensive than the computation of G 1/2 since, as in the one-dimensional case (p = 1), (21) will allow us to make use of the Fast Fourier Transform (FFT) which considerable reduces the complexity. Now, the algorithm proposed by Wood and Chan may be described through the following steps. Let m be a power of 2 greater than 2(n -1).

Step

1. For 1 ≤ u ≤ v ≤ p calculate for k = 0, . . . , m -1 B u,v (k) = m-1 j=0 C u,v (j)e -2iπjk/m
where C u,v (j) if the element (u, v) of the matrix C(j) defined by [START_REF] Lavancier | Covariance function of vector self-similar process[END_REF]

and set B v,u (k) = B u,v (k) * .
Step 2. For each j = 0, . . . , m -1 determine a unitary matrix R(j) and real numbers ξ u (j) (u = 1, . . . , p) such that B(j) = R(j) diag(ξ 1 (j), . . . , ξ p (j)) R(j) * .

Step 3. Assume that the eigenvalues ξ 1 (j), . . . , ξ p (j) are non-negative (see Remark 11) and define B(j) = R(j) diag( ξ 1 (j), . . . , ξ p (j)) R(j) * .

Step 4. For j = 0, . . . , m/2 generate independent vectors U(j), V (j) ∼ N p (0, I) and define

Z(j) = 1 √ 2m × √ 2U(j)
for j = 0, m 2 U(j) + iV (j) for j = 1, . . . , m 2 -1, let Z(mj) = Z(j) for j = m 2 + 1, . . . , m -1 and set W (j) := B(j)Z(j).

Step 5. For u = 1, . . . , p calculate for k = 0, . . . , m -1

∆X u (k) = m-1 j=0
W u (j)e -2iπjk/m and return ∆X u (k), 1 ≤ u ≤ p, k = 0, . . . , n -1 .

Step 6. For u = 1, . . . , p take the cumulative sums ∆X u to get the uth component X u of a sample path of a mfBm. Figure 3 gives some examples of sample paths of mfBm's simulated with this algorithm. Remark 11. The crucial point of the previous algorithm lies in the non-negativity of the eigenvalues ξ 1 (j), . . . , ξ p (j) for any j = 0, . . . , m -1. In the one-dimensional case (when p = 1) Steps 2 and 3 disappear, and in Step 1, B 11 (k) corresponds to the k-th eigenvalue of the circulant matrix C 11 with first line defined by C 11 (j) = γ 11 (j) for 0 ≤ j ≤ m/2 and γ(mj) for j = m/2 + 1, . . . , m -1. For the fractional Gaussian noise, it has been proved by Craigmile [START_REF] Craigmile | Simulating a class of stationary gaussian processes using the davies-harte algortihm, with application to long memory processes[END_REF] for H < 1/2, and by Dietrich and Newsam [START_REF] Dietrich | Fast and exact simulation of stationary gaussian processes through circulant embedding of the covariance matrix[END_REF] for H > 1/2 that such a matrix is semidefinite-positive for any m (and so for the first power of 2 greater than 2(n-1)). In the more general case p > 1, the problem is much more complex: the quantities B u,v (k) are not necessarily real, and the establishment of a condition of positivity for the matrix B u,v (k) does not seem obvious. When the condition in Step 3 does not hold, Wood and Chan suggest to either increase the value of m and restart Steps 1,2 or to truncate the negative eigenvalues to zero which leads to an approximate procedure. These problems are not addressed in this paper. Let us assert that for the simulation examples presented in Figure 3, we have observed that this condition is satisfied for m equal to the first power of 2.

Figure 2 represents |ρ 1 , 2 |

 12 Figure 2 represents |ρ 1,2 | with respect to (H 1 , H 2 ).Figures1 and 2illustrate the main limitation of the mfBm model. Under self-similarity condition (1), it is not possible to construct arbitrary correlated fractional Brownian motions. For example, when H 1 = 0.1 and H 2 = 0.8, the correlation cannot exceed 0.514.

Figure 1 .

 1 Figure 1. Various examples of possible values for(ρ 1,2 , η ′ 1,2 ) with η ′ 1,2 := η 1,2 (1 -H 1 -H 2 ) when H 1 + H 2 = 1 and (ρ 1,2 , η1,2 ) when H 1 + H 2 = 1,ensuring that Σ(s, t) is a covariance matrix function in the particular case p = 2.

Figure 2 .

 2 Figure 2. Maximal values of the absolute possible correlation parameter |ρ 1,2 |ensuring that Σ(s, t) is a covariance matrix function in the case p = 2, in terms of H 1 and H 2 for the causal and well-balanced mfBm.

Remark 10 .

 10 Let us discuss the computation cost of the most expensive steps, that is steps 1, 2 and 5. Step 1 requires p(p+1) 2 applications of the FFT of signals of length m, Step 2 needs m diagonalisations of p × p Hermitian matrices and Step 5 requires p applications of the FFT of signals of length m. Therefore, the total cost, κ(m, p) equals κ(m, p) = O p(p + 1) 2 m log m + O(mp 3 ) + O(pm log m).

Figure 3 .

 3 Figure 3. Examples of sample paths of length n = 1000 normalized causal (top), well-balanced (middle) and general (bottom) multivariate fractional Brownian motion with p = 5 components. The Hurst parameters are equally spaced in [0.3, 0.4] (left), [0.8, 0.9] (middle) and [0.4, 0.8] (right). The correlation parameters ρ i,j are all set to 0.4 (left, middle) and 0.3 (right). For the general mfBm (bottom), the parameters η i,j are set to 0.1/(1 -H i -H j ) and the parameters η i,j to 0.1. Note that the existence condition discussed in Proposition 9 and the condition in Step 3 of the algorithm are satisfied for these different choices of parameters. For convenience, the sample paths of the left column have been artificially shifted for visibility.

Research supported in part by ANR STARAC grant and a fellowship from Région Rhône-Alpes (France) and a Marie-Curie International Outgoing Fellowship from the European Community. Research supported in part by ANR InfoNetComaBrain grant.


Benoît Trouette

Eric Chénier

Claudine Delcarte

Béatrice Guerrier

Numerical study of convection induced by evaporation in cylindrical geometry

Numerical simulations of convection induced by solvent evaporation during the drying of a polymer solution are considered. This paper focuses on the transient thermal regime occurring at the beginning of the drying and transient solutal effects are not taken into account. The onset of convection (Bénard-Marangoni and Rayleigh-Bénard) is studied for a large range of initial thicknesses and viscosities. Several stochastic models are compared to analyze the influence of the perturbation description on the transition thresholds. Two-dimensional (2D) and three-dimensinal (3D) models are shown to give close results. The 3D model is used to characterize the pattern evolution during the drying. In the case of surface tension driven convection, a method is developed to describe the cells morphology and their time evolution.

Introduction

The drying of binary fluid films including one volatile component appears in many industrial processes such as painting, inking or packaging. When drying a polymer/solvent solution, the vaporization of the solvent induces a temperature decrease at the free surface due to the latent heat transfer and the increase of the solute concentration. The resulting thermal and mass gradients create volume (buoyancy) and surface (capillary) forces that can generate a flow into the bulk [1][2][3][4][START_REF] Colinet | Nonlinear Dynamics of Surface-Tension-Driven Instabilities[END_REF][START_REF] Merkt | [END_REF][7][8][9]. One important point is that the flow is transient in most configurations. Indeed starting from a fluid at rest at ambient temperature, a transient diffusive regime takes place as soon as evaporation begins. Over time, the heat and mass transfers evolve and can generate convection. Sometime after, the fluid motion fades away due to thermal or solutal gradients changes and/or thickness reduction and viscosity increase. At the end, transfers are governed by diffusion again.

The numerical approach presented here follows an experimental study dealing with drying of polymer solutions. The development and evolution of convective patterns induced by solvent evaporation in a Polyisobutylene/toluene solution have been experimentally studied for a large range of initial thicknesses and viscosities [10]. In this paper we focus on the thermal transient regime that concerns the beginning of the drying only; solutal effects are disregarded. 2D numerical simulations [11] and linear stability analysis [12] have been performed for the transient thermal regime. Critical Marangoni and Rayleigh numbers have been determined as a function of the experimental control parameters. As pointed out in these previous works, it is more proper to view the results as the definition of a transition region rather than well defined threshold values between stable and unstable domains. Indeed, in transient problems, the thresholds depend on the criterion used to characterize significant convection, and also on a e-mail: eric.chenier@univ-paris-est.fr the initial perturbation. It is then important to check that this transition region is thin and slightly depends on the perturbation characteristics (structure and amplitude). In the previous 2D numerical study [11], the initial condition consisted of a random vector field defined on collocation nodes. One of the purposes of this paper is to introduce more sophisticated models, closer to realistic physical behaviors, by using stochastic approaches and by taking into account spatial correlations for disturbances. The second point developed in this study is the extension of the 2D thermal model to the 3D configuration. Convection/diffusion threshold values are discussed and compared with 2D results. A particular attention is paid to flow structures at the beginning of the drying that cannot be described in 2D simulations.

Physical and mathematical model

The numerical simulations refer to experiments of polymer solutions drying [10]. Details on the approximations that have been considered for the thermal transient regime occurring at the beginning of the drying can be found in [11] and [12]. They are briefly recalled below.

A cylindrical container of radius R, with adiabatic solid surfaces, is filled with a Polyisobutylene/toluene solution over a thickness e (Fig. 1 of I.R. camera show some convective instabilities occuring on short periods of the order of some thermal diffusion times, this last time scale being 10 -3 smaller than the solutal diffusion time scale. During this thermal regime, the solutal effects and decrease of the solution thickness are neglected. The mass fraction of the solvent is then assumed constant and the viscosity used in the simulations corresponds to the initial solute concentration. The two control parameters chosen in the experiments were the initial thickness (0.3 mm ≤ e ≤ 14.3 mm) and the initial viscosity through the initial polymer mass fraction (0% ≤ ω p ≤ 15% corresponding to 5.5 × 10 -4 Pa.s ≤ µ ≤ 2.4 Pa.s). Given the physical properties of the solution and the geometric configuration, it can be shown that the instability mode due to the coupling between the Marangoni effect and the deflection of the free surface can be neglected [12]. Moreover, a non-zero base evaporation rate stabilizes the interface against deformation due to vapor dynamics [8]. The modelisation of the gas phase is here ignored and the free surface deformation is disregarded. The interface between the solution and surounding gas is then considered as flat. The fluid is assumed Newtonian. Except for the viscosity, physical properties are those of the pure solvent. They are kept constant except for both the density in the buoyant term (Boussinesq approximation) and the surface tension in the capillary stress which depend linearly on the temperature. The thermal regime is thus governed by the Navier-Stokes and energy equations in dimensionless form:

           ∂u ∂t + (u • ∇) u = -∇p + Pr∆u + RaPrΘ1 z ∂Θ ∂t + (u • ∇) Θ = ∆Θ ∇ • u = 0 (1)
where u = (u • τ )τ + (u • 1 z )1 z , p and Θ respectively denote the dimensionless velocity, pressure and temperature fields; 1 z stands for the upward unit vector normal to the free surface and τ is unitary, such that τ •1 z = 0. The position vector x writes in cylindrical coordinates (1 r , 1 φ , 1 z ), x = (r, φ, z) and in 2D-Cartesian coordinates (1 x , 1 z ), x = (x, z); the velocities are expressed by u = u1 r + v1 φ + w1 z and u = u1 x + w1 z , respectively. The length, velocity, pressure and time are scaled by the fluid thickness e, the thermal diffusion velocity α/e, the dynamical pressure ρα 2 /e 2 and the thermal diffusion time e 2 /α, where α and ρ are the thermal diffusivity and the density of the solution. The reduced temperature is defined by Θ = (T -T 0 )/∆T s where ∆T s is the difference between the initial temperature T 0 and the steady and uniform temperature at the end of the transient thermal regime. The temperature scale results from the balance, at the free surface, between the heat flux lost by evaporation, the convective heat transfer with the surrounding gas and the conduction in the liquid film. A one-sided model is used and heat transfer with the surroundings is taken into account by a constant and uniform convective coefficient h. Assuming small temperature variations, the evaporation flux can be linearized, given [11]:

-λ ∂T ∂z (z=e) = h(T (z = e) -T 0 ) + L Φ(T 0 ) + ∂Φ ∂T T0 (T (z = e) -T 0 )
with λ the thermal conductivity, Φ(T 0 ) the evaporative heat flux at T 0 , L the latent heat.

Setting

H = h + L ∂Φ ∂T | T0 , the global exchange coefficient, we obtain -λ ∂T ∂z (z=e) = H(T (z = e) -T 0 ) + LΦ(T 0 )
At the end of the transient thermal regime, the temperature becomes steady and uniform again: the heat flux gained by convective transfer with the surroundings balances exactly the heat flux lost by evaporation. The temperature difference between the liquid and the surroundings is then ∆T s = LΦ(T0)

H

. In cylindrical coordinates, the dimensionless boundary conditions write:

at r = A or z = 0    u = 0 (no slip) ∂Θ ∂n = 0 (adiabatic surfaces) and at z = 1            ∂(u • τ ) ∂z = -Ma ∂Θ ∂τ (thermocapillary stress) w = 0 (cf. [12]) - ∂Θ ∂z = Bi (Θ + 1) (energy balance) (2)
The dimensionless parameters A = R/e, Pr = ν/α, Ra = gβ∆T e 3 /µα, Ma = -∂γ ∂T e∆T /µα and Bi = H/α are the aspect ratio of the liquid film and the Prandtl, Rayleigh, Marangoni and Biot numbers, µ, ν, γ being the dynamic and kinematics viscosities and the surface tension coefficient respectively. The dimensionless parameters Pr, Ra, Ma and Bi depend on two physical control parameters used in experiments, namely the initial fluid thickness and the initial viscosity (which is a function of the polymer mass fraction [10]); the relations between control Ma Ra Pr Bi 5850 e/µ 451 e 3 /µ 12 µ 0.2 e Table 1. Relations between the dimensionless parameters and the thicknesses ([e] = mm) and viscosities ([µ] = mPa.s) parameters and dimensionless numbers are given in Table 1.

For the 3D configuration, a pseudo-spectral method has been developed to solve the equations (1) with boundary conditions (2) in cylindrical coordinates (r, φ, z) [13]. The solution in the azimuthal direction has been expanded in Fourier series with 128 complex conjugate modes. The evaluations of the radial and axial space derivatives are based on Chebyshev collocation method with typically 181 Gauss-Radau and 20 Gauss-Lobatto points. Convective and advective terms are handled with a second order Adams Bashforth scheme, the diffusion contributions being implicitly treated. Velocity and pressure decoupling is achieved by projection algorithm. Computations are initialized with the fluid at rest and with small disturbances added to the uniform temperature field (cf. next section).

To study the sensitivity of the transient flow to the initial disturbances, a 2D model has also been used. The equations have been solved, in Cartesian geometry (x, z), using finite volume discretization with second order space approximations. The rectangular liquid layer length 2A corresponds to the diameter of the container. The time discretization scheme and the pressurevelocity decoupling are similar to those used in 3D simulations.

Uncertainty model for the initial thermal condition

The stochastic approach proposed in this paper aims to introduce relevant uncertainties in the model in order to better take into account experimental conditions.

As previously said, the fluid is initially at rest, u = 0, and perturbations are applied on the initial temperature only. The mean value of the resulting field is equal to the ambient temperature, i.e. zero in reduced variables. Using the subscript ξ to denote stochastic variables, the initial temperature field Θ ξ (x, t = 0) writes:

Θ ξ (x, t = 0) = 0 + θ ξ (x)
The random variables ξ are here provided by a uniform and centered probability density function p(ξ) = √ 3/6 on [-√ 3;

√ 3] such that the mean value E(ξ) = R ξp(ξ) dξ = 0 and the variance E(ξ 2 ) = R ξ 2 p(ξ) dξ = 1. The evolutions of the velocity and temperature are then analyzed by means of mean values and standard errors. These statistical quantities are computed by a Monte Carlo method which necessitates numerous deterministic simulations. The random samples are obtained by a Twister Mersenne algorithm [14]. Several stochastic models of θ ξ (x) have been compared. The first model hes been used for both 2D and 3D simulations and the others only in the 2D approach in Cartesian geometries, in order to limit computational time. For simplicity, equations are only given for the 2D-Cartesian configuration, with -A ≤ x ≤ A and 0 ≤ z ≤ 1.

1. The simplest idea consists in taking into account the uncertainties directly in the numerical scheme. To do this, we define the discrete initial condition so that a random value, defined by the probability law p(ξ), is accociated to each collocation node for the temperature. Using the finite volume method with

N x +1 discrete x-coordinates x i = (2i-N x )/N x ×A, i = 0, • • • , N x and N z +1 discrete z-coordinates z j = j/N z , j = 0, • • • , N z ,
and two independent families of random variables {ξ i , i = 0, • • • , N x } and {ξ j , j = 0, • • • , N z } defined by p(ξ), the stochastic initial condition writes:

θ ξ (x, z) = a √ 3 6 Nx-1 i=0 Nz-1 j=0 H(x -x i )H(x i+1 -x)H(z -z j )H(z j+1 -z)ξ i ξ j (3)
with H(x) the Heaviside function and a the maximum magnitude of the disturbance. Such an approach is very easy to implement but may be questionable because of the non regularity of the initial thermal condition. Indeed spectral methods, and in a smaller extent finite volume methods, require regular solutions. The gradients of the initial temperature diverge as a function of the mesh size, so that convergence analysis cannot be easily performed. Despite this point, this simple stochastic approach has been used for all computations performed with the spectral 3D code. 2. In order to guarantee the regularity of the initial thermal condition, we assume that disturbances are spatially correlated. This assumption is also physically relevant because the smallest patterns are necessarily driven by diffusion scales. For the sake of simplicity and considering the lack of information about physical disturbances, we suppose that the random thermal initial condition only depends on the horizontal coordinate x and is invariant in the z-direction. The stochastic initial condition is then expressed by a Karhunen-Loève decomposition [START_REF] Rozanov | Processus aléatoires[END_REF]:

θ ξ (x, z) = ∞ m=0 β m ξ m φ m (x) (4)
where {φ m (x)} is a Hilbert base in L 2 ([-A; A]), β m > 0 such that β i+1 ≤ β i and lim n→∞ β n = 0, and {ξ m } is a family of random values defined by the probability law p(ξ). The autocorrelation function R(x, x ′ ) was chosen to enable analytical calculations:

R(x, x ′ ) = R 0 exp -|x -x ′ | λ
with λ the correlation length [START_REF] Chénier | Sensitivity of diffusive-convective transition to the initial conditions in a transient Bénard-Marangoni problem[END_REF]. It can be shown that:

                         φ n (x) = A n cos(ω n x + χ n ) A n = 2ω n 2Aω n + sin(2Aω n ) cos(2Aω n + 2χ n ) β n = 2λR 0 1 + λ 2 ω 2 n χ n = -Aω n + nπ 2 tan Aω n - nπ 2 = 1 ω n λ (5)
3. The third stochastic model uses an uncorrelated disturbance similar to the first model, but only in x-direction. Like in the second model the random thermal initial condition is therefore invariant in z-direction.

θ ξ (x, z) = a √ 3 6 Nx-1 i=0 H(x -x i )H(x i+1 -x)ξ i (6)
In order to introduce the same thermal energy (L 2 -norm) at the initial time, the intensity of the auto-correlation function fulfills the equality R 0 = a 2 /12.

Results

Onset of convection

The results presented in this paragraph are obtained for e = 1 mm and µ = 3 mPa.s and with A = 10. Statistics are achieved with 100 and 200 random shots in 3D and 2D, respectively. Fig. 2 shows the transient behavior of the fluid and the comparison between the different stochastic models. Filled symbols and crosses stand for 2D-and 3D-simulations, respectively; 3)), triangles (Eqs. ( 4) and ( 5) with λ = 1) and squares (Eq. ( 6)) are 2D-results. For each symbol set, the area between the two curves corresponds to the uncertainty domain based on the mean value, plus and minus the standard error.

curves with the same symbols correspond to mean value plus or minus standard error. Fig. 2(a) presents ∆Θ max , the maxima of the temperature difference between the bottom and top of the fluid layer as a function of time. It must be emphasized that, when convection is on, these maxima are larger than in the pure diffusive regime (full line curve). This is linked to local effects. Indeed, once the flow is convective, the vertical thermal gradient at the free surface decreases in the regions where the fluid dives into the core bulk. The conduction being reduced at the interface (cf. Eq. ( 2)), the temperature at the free surface decreases due to the heat lost by evaporation (cf. negative values at z = 1 in Fig. 3). In contrast at the bottom, the temperature is larger than in the diffusion case when the fluid flows upward (cf. positive values at z = 0 in Fig. 3). Consequently, the maximum temperature difference between the bottom and the interface increases beyond the diffusive value. This effect is more pronounced in the 3D model since the Peclet number is about 40% higher than in the 2D case. However the mean value of the temperature difference is always smaller in the convective regime than in the diffusive regime, as expected. The curves depart from the diffusive case when the convection becomes significant ( 0.1 < t < 0.3). These times correspond to a Peclet number close to one, i.e. when the thermal diffusion velocity and the maximum of the velocity are of the same order (Fig. 2(b)). Several conclusions can be drawn from these results. Whereas the mean initial thermal energy was the same in the different stochastic models, convection arises earlier with correlated disturbances because the correlation length, λ = 1, is close to the natural wavelength of the convective cells (cf. Fig. 3). Disturbances based on uncorrelated random values (Eq. ( 3)) are the less effective. However the use of the three stochastic models does not affect significantly the transient flow and heat transfer, what is in good agreement with [START_REF] Chénier | Sensitivity of diffusive-convective transition to the initial conditions in a transient Bénard-Marangoni problem[END_REF] and justifies the use of the simple first stochastic approach. 2D and 3D simulations give rather close results. Although the maximum value of the Peclet number is 40% higher in the 3D case, the departure from the diffusive regime occurs more or less at the same time.

In a second step the influence of the perturbation amplitude on the onset of convention is studied. To characterize the presence of observable convection, a criterion based on the thermal Peclet number is chosen. Convection will be considered as significant if, when the system is submitted to an initial perturbation, there is a time t at which the perturbation is significantly amplified, i.e. such as Pe = max x |u(x, t)| > 1. Results are presented in the viscosity-thickness map (Fig. 4), for the 3D model and for perturbation amplitudes varying on several orders of magnitude,

√ 3 • 10 -9 ≤ a ≤ √ 3 • 10 -3 (cf.
Eq. ( 3)). As a comparison the thresholds computed with the 2D model by Touazi et al. [11] are also displayed. As is made clear in the figure, on 3D-results 2D-results [11] Thickness [mm] Viscosity [mP a

• s] the one hand the dispersion is small for the large explored range of perturbation amplitudes, on the other hand 2D and 3D simulations give close results.

In conclusion, these results complete and confirm previous analyses concerning the onset of convection for the thermal transient regime [11,12]. Stability analysis is meaningful for this regime, in the sense that it is possible to characterize a thin transition region in the control parameters space between stable and unstable configurations. This region is not very sensitive to disturbance description (shape and amplitude) and is similar in 2D and 3D models.

At last, some simulations have been performed on 3D configurations of smaller aspect ratio. Thresholds are not very different, since a relative difference of less than 11% is obtained when reducing the aspect ratio from A = 10 to A = 3 and, of less than 43%, for a reduction from A = 10 to A = 1.

Convective patterns

Unlike 2D simulations, the 3D configuration allows to describe more accurately the convective patterns and their time evolutions. Fig. 5 and Fig. 6 show the temperature field evolutions at the free surface for two configurations, e = 1 mm, µ = 3 mPa.s and e = 30 mm, µ = 3500 mPa.s respectively. Flow patterns are different for these two configurations. Indeed for small thickness corresponding to surface tension driven flow, cells are observed as in classical Rayleigh-Bénard-Marangoni studies [START_REF] Médale | [END_REF]. They grow in time but become also less active, since P e and ∆Θ max decrease for t > 1 (cf. Fig. 2 and Fig. 5). Additional patterns emerge for the largest thickness linked to buoyancy driven flow. The convective flow organizes itself into rolls (Fig. 6) [10]. Fig. 6. Iso-temperaturelevels at the interface, for several successive times, e = 30 mm, µ = 3500 mPa.s. The Peclet number Pe and maxima of temperature differences in the fluid layer ∆Θmax and at the free surface ∆Θ fs are specified.

To get a more systematic characterization of capillary cells for e = 1 mm and µ = 3 mPa.s, a reconstruction method has been developed. First, the positions the temperature maxima on the free surface are detected. Then the cell faces are defined as belonging to perpendicular bisectors between the nearest adjacent hot points. This analysis is performed with a Voronoi algorithm (cf. Fig. 5). Fig. 7 presents the evolution of the mean value of the cells size with its standard error. The mean cell size increases by a factor of about two while, as said above, convection fades away (Fig. 2(b)). For purpose of comparison, 2D results, for the same mean initial thermal energy, are also reported. The mean 2D wavelengths, λ 2D , are defined by the couples of contra-rotating convective cells (λ 2D = (2A)/(2n) with 2A the cavity extension and n the number of cells). For t > 0.5, the mean sizes of 3D-capillary cells at the free surface are in good agreement with the 2D patterns. The 2D model seems to overestimate the wavelength of about 30% at the beginning of the drying. However, it is important to recall that the flow is then almost conductive and it is thus very difficult to characterize the cell size.

To go a little further in the characterization of the cell shapes and their time evolution during the transient thermal regime, cells coming from 100 different shots are classified with respect to the number of edges as follows. For each cell, only edges greater than 0.1 times the largest side are considered. Cells in contact with the boundary are not taken into account. It can be seen in Fig. 8 that the majority of the patterns are pentagons and hexagons and the proportions do not evolve noticeably. The hexagonal shape is very usual in the classical Bénard-Marangoni configuration, provided that the aspect ratio of the liquid layer is large enough and the cells are far from lateral walls. Notice that pentagons and squares have also been observed in numerical and experimental studies in case of square containers [START_REF] Médale | [END_REF]. The present categorization of flow structures is only based on geometrical considerations but the observed patterns can differ substantially from their ideal shapes (edges of same length) as obtained in the steady Bénard-Marangoni problem [START_REF] Médale | [END_REF]. To measure the departure of each cell from its regular polygon shape, the following procedure is used. For a hexagonal cell for instance, we measure its area and define the reference regular hexagon having the same area. The cell is divided up into six triangles with the cell center as apex and the six cell edges as base of the triangles. Then the standard error σ between the area of the reference triangle (1/6 of the hexagon area) and the six triangles of the cell is computed and used to characterize the cell regularity. Results are given in Fig. 9. The percentage of cells having a standard error in the range [20%; 40%] is about 60% all along the evolution of the transient regime. The cells induced by evaporation with random initial conditions are less regular than those observed in the usual steady Bénard-Marangoni configuration. Notice however that, more regular cells are recovered if initial conditions are smooth or well structured in space, but this necessitates a very accurate control of the initial condition, what is really hard to achieve in experimental setups.

Conclusion

The thermal transient regime, occurring at the beginning of the drying of a polymer/solvent solution, has been studied numerically. A particular attention has been payd to the impact of the initial conditions by means of several stochastic models of perturbations. The main result is that it is possible to characterize a thin transition region between stable and unstable configurations and this region is not very sensitive to the shape and amplitude of the disturbance. The onset of convection has been determined, with 2D and 3D models, for a large range of initial thicknesses and viscosities. The results complete and confirm previous 2D stability analysis. 2D and 3D thresholds are not very different. The 3D model has also been used to characterize the flow pattern evolution during the drying. In the case of convection driven by surface tension, tools have been developed in order to accurately describe the cell morphology. It is shown that cells are mostly non-regular pentagons or hexagons. The cell size increases during the drying while convection fades away. The next step of this study will concern the transient solutal convection. This problem is more complex since it implies to take into account both thickness reduction and viscosity increase.

Fig. 1 .

 1 Fig. 1. Scheme of the configuration

Fig. 2 .

 2 Fig.2. Transient evolutions for e = 1 mm, µ = 3 mPa.s and A = 10. The full line corresponds to the diffusion case (Ra = Ma = 0). Crosses stand for 3D-results with A = 10. Filled circles (Eq. (3)), triangles (Eqs. (4) and (5) with λ = 1) and squares (Eq. (6)) are 2D-results. For each symbol set, the area between the two curves corresponds to the uncertainty domain based on the mean value, plus and minus the standard error.

Fig. 3 .

 3 Fig. 3. Difference between the convective and conductive thermal fields and equally spaced streamlines, for e = 1 mm, µ = 3 mPa.s, -10 ≤ x ≤ 0 and t = 1

Fig. 4 .

 4 Fig. 4. Diffusion/convection transition as a function of the liquid thickness and viscosity, A = 10. In the 3D-model, the perturbations amplitude a (cf. Eq. (3)) ranges from √ 3 • 10 -9 to √ 3 • 10 -3 .

  (a) t = 0.7, Pe = 22.95, ∆Θmax = 12 × 10 -2 , ∆Θ fs = 7.44 × 10 -2 (b) t = 1, Pe = 20.33, ∆Θmax = 13.2 × 10 -2 , ∆Θ fs = 8.73 × 10 -2 (c) t = 6, Pe = 7.77, ∆Θmax = 6.26 × 10 -2 , ∆Θ fs = 3.55 × 10 -2 (d) t = 8, Pe = 3.47, ∆Θmax = 3.57 × 10 -2 , ∆Θ fs = 1.56 × 10 -2

Fig. 5 .

 5 Fig.5. Iso-temperature levels and reconstructed cells at the interface, for several successive times, e = 1 mm, µ = 3 mPa.s. The Peclet number Pe (cf. Fig.2(b)), the maxima of temperature differences in the fluid layer ∆Θmax (cf. Fig.2(a)) and at the free surface ∆Θ fs are specified.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Dimensionless cell size: mean values and standard errors

Fig. 9 .

 9 Fig. 9. Evolution of the regularity measured by standard error with respect to regular pentagons or hexagons. e = 1 mm, µ = 3 mPa.s.

Acknowledgment

This work has been supported by the "Fédération TMC (Ministère de l'Enseignement supérieur et de la Recherche, France)"