Nishanth Dongari
email: nishanth.dongari@strath.ac.uk

Yonghao Zhang

Jason M Reese

Molecular free path distribution in rarefied gases

Keywords: numbers: 47.45.-n, 47.61.Cb, 47.11.Mn

We present the results of investigations into the distribution of molecular free paths in rarefied gases using molecular dynamics simulations. Our tests on a range of different gas densities and confinements (unbounded, single bounding wall, and parallel bounding walls) indicate that the molecules perform Lévy-type flights, irrespective of the presence of a bounding wall. The free paths most closely follow a power-law distribution. Simulations of gases confined by planar surfaces indicate that the local molecular mean free path varies sharply close to a solid surface. These results may yield new insight into diffusive transport in rarefied gases, in particular, the constitutive behaviour of gas flows in micro and nanoscale devices.

Introduction

Rarefied gases often display non-standard fluid behaviour, as demonstrated in numerous experiments [START_REF] Karniadakis | Microflows and Nanoflows: Fundamentals and Simulation[END_REF]. The key parameter indicating the rarefaction of the gas (and the degree of departure from quasi-equilibrium fluid mechanics) is the Knudsen number, Kn, defined as the ratio of the gas molecular mean free path λ to the characteristic system length scale L. Understanding the behaviour of gases over a wide range of Kn is key to the design of future micro and nanoscale devices and to the manipulation of materials at the molecular scale.

While bulk gas flow properties, such as the average mass flowrate along a micro channel, can be experimentally measured, interpretation of a single integral flow parameter in terms of the components making up the flow model is problematic. Any deviation of an integral quantity from its value calculated using conventional macroscopic fluid mechanics could be ascribed to multiple and competing causes, e.g. deficiencies in the bounding surface slip model, the effect of Knudsen layers, non-linear constitutive behaviour in the bulk, etc. A posteriori fitting of macroscopic theories to the experimental data with the aid of tuning parameters is unlikely to be physically realistic, or general enough, although it is widely practised [START_REF] Karniadakis | Microflows and Nanoflows: Fundamentals and Simulation[END_REF].

In this situation, and in order to shed light on the molecular phenomena that determine the macroscopic behaviour, simulating rarefied gases determinstically comes to the fore. Simulation methods such as Molecular Dynamics (MD) can allow us to assess fundamental molecular properties of the gas -the inter-molecular collision rate, the velocity distribution, the distribution of free paths between collisions -over a range of rarefied conditions, and in both unconfined and confined spaces [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF]. MD is the most appropriate method as it is deterministic, allowing for realistic molecular behaviour, i.e. molecular attractions, repulsions, movements and scatterings.

This paper presents the results of numerical MD experiments on simple rarefied gases. We highlight the substantively different behaviour of the molecular free paths, collision rates, and mean free path profiles from their classical equilibrium counterparts. In particular, we show that Lévy-type or power-law descriptions better describe the rarefied gas behaviour. As the diffusive transport coefficients of gases can be interpreted in terms of the collisions of molecules, and of the free paths of the molecules between collisions [START_REF] Kennard | Kinetic Theory of Gases with an Introduction to Statistical Mechanics[END_REF][START_REF] Cercignani | Rarefied Gas Dynamics: from Basic Concepts to Actual Calculations[END_REF], our results may yield new insight into the fluid mechanics of the rarefied case. This work may also be relevant to the construction of particle-mesh simulations (e.g. the DSMC method) of rarefied flows in micro and nano systems.

Methodology

Molecular dynamics modelling

We use the open source software OpenFOAM [START_REF] Openfoam | The Open Source CFD Toolbox[END_REF], which includes the MD routines implemented by Macpherson et al [START_REF] Macpherson | [END_REF]7], although modified as we describe below. Monatomic Lennard-Jones molecules are simulated, such that the short-range repulsive/long-range attractive potential between two molecules separated by a distance r is:

V LJ (r) = 4 r σ -12 - r σ -6 , (1)
where is related to the interaction strength of the molecules and σ is the hardsphere diameter of the gas molecules [START_REF] Cercignani | Rarefied Gas Dynamics: from Basic Concepts to Actual Calculations[END_REF] (see also section 2.2 below). We investigate Neon gas, which has σ = 2.5830 × 10 -10 m, = 7.1957 × 10 -22 J, and a molar mass m = 20.180 kg/kmol. The molecules are spatially distributed in a domain of interest with a random initial Gaussian velocity distribution corresponding to a prescribed gas temperature [START_REF] Kennard | Kinetic Theory of Gases with an Introduction to Statistical Mechanics[END_REF]. They are then allowed to relax to their steady-state configuration before we take measurements.

All interactions between nonbonded particles should be calculated but, since the LJ potential vanishes at larger r, only the interactions between molecules within a certain cutoff radius r c need to be evaluated in the MD simulation. We choose a cut-off radius of r c = 2.5σ, as any larger value has a negligible effect on our results but increases the computational cost of the simulations. The equations of molecular motion are integrated using a leapfrog scheme [7] with a step size of 0.005τ , where τ = (mσ 2 / ) (∼ 10 -12 s). The spatially-resolved flow parameters are written at intervals of 20τ , and the total run time of the simulations to relax to steady-state is at least 3000τ (in the extreme rarefied case we test, up to 50000τ ).

In order to record the free paths of individual molecules in MD simulations, conventionally each molecule's travelled distance between two successive collisions is required. But this can only be recorded when collisions occur and, in simulations with periodic boundaries, it is not appropriate to directly measure the travelled distance of the molecule as a difference between its current position and the position of its last collision. Instead, we obtain the molecular free path as l = (t C -t LC )Sp C , where t C is the current simulated time, t LC is set to the current simulation time when a collision occurs (so it holds the time of the last occurring collision), and Sp C is the current speed of the molecule. These molecular attributes, t LC and Sp C , we implemented into the MD solver in OpenFOAM, along with the attribute col that keeps track of whether a target molecule is currently undergoing a collision. The latter is activated (set to 1) when the target molecule is within r col of any other molecule (see section 2.2 below). It is deactivated (set to 0) if it is currently activated and the target molecule is greater than r col from any other molecule.

This new method records the molecular free paths in every time step, and its major advantage is that we can calculate the probability that molecules travel a distance without experiencing a collision, i.e. the free path distribution function. The simple average of these free paths of all simulated molecules is the mean free path. Additionally, the collision handling algorithm allows the evaluation of multiple collisions, i.e. collision events which involve three or more molecules simultaneously. This means the dilute gas assumption employed in deriving the Boltzmann equation can be explored.

The collision radius

We need to identify the occurrence of inter-molecular collisions in the MD simulations. While the instance of a collision is unambiguous only if the molecular interaction is a step potential (i.e. hard-sphere molecules), the time step required for simulations in that case would need to be vanishingly small in order to capture the collision instant reliably [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF]. For the continuous and more physically-realistic LJ potential that we use, we instead declare a collision to have occurred if two molecules are closer together than a distance r col . Molecules in the simulation start to exchange momentum as soon as they are within the MD cut-off distance r c , and they start to repel each other when they are closer than σ. Hence, the closeness parameter r col must lie between 2.5σ and σ.

An appropriate value for r col can be evaluated through MD simulations to obtain the mean free path for unconfined Neon gas -λ Ne = (µ/ρ)(π/2RT ) 1/2 , where µ is the gas dynamic viscosity, ρ is gas density, R is specific gas constant and T is temperature -at standard (i.e. non-rarefied) conditions. While computationally time-consuming at standard gas densities, this is also an effective way of validating the MD technique overall. At STP, the expression above give a λ Ne of 1.233 × 10 -7 m [START_REF] Cercignani | Rarefied Gas Dynamics: from Basic Concepts to Actual Calculations[END_REF]. Using r col = σ, the MD simulation retrieves this mean free path value within 1% (see section 3.1.1). As a sensitivity check, if we use r col = 2.5σ our MD simulation result is 15% less than the theoretical mean free path value. So the appropriate value for r col is σ, and all the simulations in this paper have been carried out using this value.

Results and Discussion

Unconfined gases

Collision rate probabilities, the distribution of free path lengths, and the steady state mean free path (MFP) for different rarefied gas conditions are simulated using MD in an unconfined cubic domain.

Mean free path

To validate our MD method, and as noted in section 2.2 above, we first use it to deduce the unconfined value of the mean free path of Neon gas at standard conditions. We simulate a cubic spatial geometry with periodic boundary conditions, i.e. if a molecule leaves the domain, it is reintroduced on the opposite side of the domain, retaining the same velocity. The side-length of this simulation domain is chosen to be the MFP of Neon at STP conditions, and the domain is filled with 51505 molecules to achieve a gas density of 0.9011 kg/m 3 .

Figure 1(a) shows the recorded MFP as the simulation progresses in time. The MD measurements relax to an MFP value of 1.221× 10 -7 m, that differs from the theoretical value of MFP by only 0.96%. The simulations also recorded the number of binary intermolecular collisions (by discounting a molecule as undergoing a binary collision if its counterpart molecule is already undergoing a collision, i.e. a molecule with col = 0 cannot collide with a molecule with col = 1). At standard density, binary collisions account for 96.92% of all collisions, the remainder being multiple-molecule collisions. If only binary collisions are allowed to contribute to the MFP measurement we retrieve an MFP of 1.262 × 10 -7 m, which is close to the kinetic value of 1.255 × 10 -7 m given through [START_REF] Kennard | Kinetic Theory of Gases with an Introduction to Statistical Mechanics[END_REF]:

λ k = m/N Av πρσ 2 √ 2 , (2)
where N Av = 6.0221415 × 10 23 is Avogadro's number. This indicates that the kinetic theory description of MFP is accurate if collisions are purely binary, as expected.

We then test a rarefied case: Neon gas of density 0.1055kg/m 3 . The side-lengths of the simulation cubic domain were increased to 2.5λ Ne , and the simulations run up to 10 × 10 -9 s of problem time because, with a reduced frequency of inter-molecular collisions, the steady-state MFP is reached more slowly. Figure 1(b) shows the MFP measurements from the simulations. As expected in this rarefied case, the equilibrium kinetic value of MFP given by (2) differs more markedly from the deterministicallysimulated MFP, in this case by about 2.31%. More than 99% of the inter-molecular collisions are binary, and fluctuations in the steady-state MFP value increased.

Collision rates and free path distributions

Besides the value of the mean free path, the distribution of individual molecular free paths is also of interest [START_REF] Kennard | Kinetic Theory of Gases with an Introduction to Statistical Mechanics[END_REF]. For a group of similar gas molecules moving at speed v with a collision rate of θv , the molecular mean free path is λ = v/ θv . If the gas molecules are in thermodynamic equilibrium, then the distribution of free path lengths, ψ(r), can be derived as [START_REF] Kennard | Kinetic Theory of Gases with an Introduction to Statistical Mechanics[END_REF]:

ψ(r) = λ -1 p(r) with p(r) = exp(-r/λ), (3)
where p(r) is the probability of a molecule traveling distance r without experiencing any collision with other moving molecules in the group. If the gas is not bounded, the MFP of the gas molecule ensemble is then simply 0 ∞ rψ(r)dr = λ. By thermodynamic equilibrium we mean that the probability of a certain microscopic state, averaged over the details of the interactions, does not change in time or space [START_REF] Cercignani | Mathematical Methods in Kinetic Theory[END_REF]. So (3) is only valid under equilibrium conditions, i.e. if θv is constant or fluctuations in θv are negligible.

Figure 2 shows the variation of collision rate per unit number density with MD simulation time, for gas densities of 0.9011 kg/m 3 and 0.1055 kg/m 3 . At atmospheric density conditions, Fig. 2(a), the mean and standard deviation of the collision rate are 4.1 × 10 -11 and 9.49 × 10 -12 , respectively, so the temporal fluctuations are low. For the rarefied case, Fig. 2(b), the fluctuations are significantly higher and the ratio of standard deviation to the mean of collision rate is around 0.77. This indicates that the thermodynamic equilibrium assumption of a collision rate constant in time is not valid in rarefied conditions. Consequently, the probability distribution of free paths is not necessarily exponential in form in the rarefied case.

Our MD results allow us to test other forms of distribution function for the rarefied case. We note that non-equilibrium transport is often better described using Lévy or power-law statistics [START_REF] Fa | [END_REF]10]. Montroll & Scher [10] pointed out that a finite moment of the probability distribution function implies an exponential character of the randomness. A distribution function with diverging higher-order moments, such as the standard deviation, is essential to non-equilibrium transport. In transient photoconductivity experiments, and the adsorption of ions in semiconductors, concentration distributions may have a long tail such that the standard deviation is diverging [11]. This longtail problem is common to various flows in non-uniform media, in many fields such as chemical engineering and environmental sciences [12].

Montroll & Scher [10] and Schlesinger [13] showed theoretically that anomalous transport is characteristic of a distribution function with diverging higher-order moments, e.g. a power-law (PL) form for the distribution function. Second moments (i.e. the standard deviation) diverge for many PL functions used to describe the probability density distribution of electrons in photoconductivity experiments.

We hypothesize that the molecular free paths in a rarefied gas follow a power-law (PL) probability distribution function. We therefore investigate the following molecular free path distribution function [16]:

ψ(r) = C(a + r) -n , (4)
where a and C are constants with positive values determined through the zeroth and first moments. A gas has a finite MFP value, so the first moment of the distribution function is required to be finite. The range of values for the exponent n can be obtained by making one of the higher-order moments divergent. Zeroth and first moments are given as follows:

1 = 0 ∞ C(a + r) -n dr, (5)
λ = 0 ∞ Cr(a + r) -n dr. (6)
Equation ( 5) requires the probability to range only from zero to one; (6) defines the unconfined MFP value. It then follows that C = (n -1)a n-1 and a = λ(n -2). If n > 2 then a > 0, otherwise the distribution function is negative. The standard deviation (second moment) of the distribution function (4) diverges only for n ≤ 3; so 2 < n ≤ 3.

If we wish to make the i th moment diverging, then max(n) = i + 1. If n → ∞, the distribution function will have finite moments, which is the condition required of an equilibrium distribution function of exponential form. For a finite n, the distribution function describes a system deviating from equilibrium. So n acts as a decisive parameter to define the extent of deviation from equilibrium. Figure 3 presents our MD results for the probability a molecule travels a distance r/λ without experiencing a collision, at both atmospheric and rarefied conditions (where λ refers to the theoretical MFP at the corresponding densities). At atmospheric density, both MD measurements and the power-law model agree well with the classical exponential distribution function. In this case n = 16 is used in the power-law model to fit the MD data, and this larger value of n conveys that the system is quite close to equilibrium. Minor deviations from the exponential distribution are noticed, and this is probably due to the consideration of multiple particle collisions.

However, in the low density case, Fig. 3(b), the deviations between the MD measurements and the classical distribution function are significant, with the discrepancies relatively higher at the longer free paths (see inset of this figure). The power-law model with an exponent n value of 9 gives good agreement with the MD data, and this decrease in n illustrates that the system is relatively more non-equilibrium.

Both the MD data and the power-law function predict a long tail to the distribution function, while classically there is a much faster decay of free paths. A faster decay indicates the system does not have many longer flights of molecules, and it follows the Brownian motion described by Einstein [START_REF] Einstein | Investigations on the Theory of Brownian Movement[END_REF]. However, this is only valid for times much larger than a characteristic time that, in turn, is much larger than the times characterising molecular collisions: this condition may not be fulfilled under rarefied conditions. The MD data show longer Lévy-type free flights of molecules, and the long tail of the distribution function implies higher-order moments are diverging.

Confined gases

The mean free path of the gas molecule ensemble approaches the theoretically predicted value only if the gas is unbounded (i.e. r = 0 → ∞). If a solid bounding surface is included in the system, some molecules will hit the surface and their free flight paths will be terminated. The MFP of all the gas molecules in the system will therefore be reduced due to this boundary limiting effect. Stops [START_REF] Stops | [END_REF] derived such a geometry-dependent MFP but used the classical exponential free path distribution function; Dongari et al [16] have developed this further using the PL distribution [START_REF] Cercignani | Rarefied Gas Dynamics: from Basic Concepts to Actual Calculations[END_REF].

We now make MD measurements of the MFP of gas molecules when bounding planar surfaces are part of the system, and compare our simulation results with other geometry-dependent MFP models [START_REF] Stops | [END_REF]16]. All our simulations are in a six-sided configuration, with two sides reflective bounding surfaces, and the other two pairs of sides periodic boundaries. This replicates a parallel plate configuration, with infinite width and length. The reflective surfaces are specular, i.e. the tangential velocity of molecules colliding with them is maintained, while the molecular normal velocity changes sign. Molecular reflections on the planar surfaces are recorded by setting t LC to the current simulated time, as for gas inter-molecular collisions.

We investigate gases confined by both single and parallel planar surfaces: a gas molecule's collision with a surface terminates that molecule's free path, so the local MFP should be shortened in the near-wall region. In the single-surface case, we use domain side lengths of λ Ne in the directions of the periodic boundaries, and 2λ Ne in the bounding surface directions. One of the reflective wall surfaces is required to simulate the bulk of the gas, so all reflected molecules at this surface are set to have experienced a collisionless travel of one mean free path. This is done by setting t LC = t C -λ Ne /Sp C . The gas is represented by 102950 simulation particles, and the sampling of the local MFP, λ eff , variation across the system is made over 2 nanoseconds of problem time, taken after 3 nanoseconds relaxation.

Figure 4 shows the variation of the normalized local MFP β = λ eff /λ Ne with normalized distance from the bounding surface y/λ Ne . The MD data agrees with the PL model for the free path (with an exponent n = 3), although minor deviations are noticed in the bulk region. Both the MD measurements and the PL results have a sharp gradient close to the surface, while the conventional exponential model has shallower gradients and underpredicts the MFP values in the wall vicinity. All three solutions agree on the MFP value in the bulk region, as expected.

For the parallel-surface case, a similar geometry is used and the Knudsen number is tuned by changing the gas density value. Results for the normalized local MFP between the two parallel surfaces are presented in Fig. 5. The total simulation time varied from 5 to 20 nanoseconds, based on the degree of rarefaction. At Kn = 0.2, the PL model with an exponent n = 3 is in fair agreement with the MD data in the near-wall region but deviates slightly in the bulk region. The exponential model underpredicts the MD data near-wall but shows fair agreement in the bulk. As the value of Kn increases the exponential model fails to predict the local MFP in the wall region, as well as in the bulk (so the exponential model results are not shown for comparison at Kn = 2). The MFP values predicted by the PL model compare very well to the MD data for both Kn = 0.5 and 1, although it overpredicts for Kn = 1 in the bulk. By Kn = 2 the PL model shows significant deviations from the simulation data and overpredicts in the near-wall region, although there is fair agreement in the bulk region. The MD data shows a relatively sharper gradient of effective MFP in the near-wall region, than the theoretical predictions.

Conclusion

The results of numerical molecular dynamics (MD) experiments support the argument that the molecular behaviour of rarefied gases is better described by Lévy or power- Comparison of molecular dynamics simulation data with power-law (PL) [16] and exponential [START_REF] Stops | [END_REF] effective local MFP models for Knudsen numbers ranging from the slip to the continuum-transition regimes. Symbols and lines are as in Fig. 4.

law types of distribution functions than the classical exponential function (or variations thereof). Inter-molecular collision rate fluctuations are significantly high in the rarefied case, which highlights the non-equilibrium nature of the gas. The MD-measured probability distribution of molecular free paths is distinguished by having a long tail that is not present in the classical distribution function.

The MD simulations in this paper are for simple geometries with planar surfaces, and for isothermal gases with no flow. It is important to explore further any new description of the molecular behaviour of rarefied gases through simulations of gases in complex geometries, with specular-diffusive and explicit walls, and also for nonisothermal cases.

Confirmation of this power-law behaviour of rarefied gases may help in the development of both new constitutive models for rarefied gases in micro and nano systems, and particle-mesh methods (such as the DSMC technique) that require a mesh that relates to the (local) molecular mean free path.

Figure 1 .

 1 Figure 1. Relaxation to steady-state of the measured mean free path as the simulation progresses, for gas densities of (a) 0.9011 kg/m 3 and (b) 0.1055 kg/m 3 . The theoretical values of the mean free path in each case are given by the straight dotted lines [4].

3 Figure 2 .

 32 Figure 2. Collision rate per unit number density as a function of simulation time, for gas densities (a) 0.9011 kg/m 3 and (b) 0.1055 kg/m 3 . The mean values of all samples are indicated by the straight lines.

5 Figure 3 .

 53 Figure 3. Probability molecules experience a collisionless travelling distance r/λ, in Neon gas of density (a) 0.9011 kg/m 3 and (b) 0.1055 kg/m 3 . Molecular dynamics measurements are compared with the exponential and power-law distribution functions. Power-law exponent n values are 16 and 9 in cases (a) and (b), respectively.

Figure 4 .

 4 Figure 4. Variation of normalized mean free path, β = λ eff /λ Ne , with normalized distance from the surface, for the planar single-surface case. Comparison of MD simulation data with power-law (PL)[16] and exponential[START_REF] Stops | [END_REF] MFP models.

2 Kn = 0. 5 Kn = 1 Kn = 2 Figure 5 .

 25125 Figure5. Variation of normalized mean free path, β = λ eff /λ Ne , with normalized distance from a surface, for the planar parallel-surfaces case. Comparison of molecular dynamics simulation data with power-law (PL)[16] and exponential[START_REF] Stops | [END_REF] effective local MFP models for Knudsen numbers ranging from the slip to the continuum-transition regimes. Symbols and lines are as in Fig.4.

Figure

  Figure 1b (DongariZhangReese_JPhysD_paper/MFP_0.1055.eps)

Acknowledgments

The authors would like to thank Matthew Borg of the University of Strathclyde, and Erik Arlemark of the Technische Universiteit Eindhoven, for very useful discussions. The research leading to these results has received funding from the European Community's Seventh Framework Programme FP7/2007-2013 under grant agreement ITN GASMEMS no. 215504. JMR would like to thank the Royal Academy of Engineering and the Leverhulme Trust for support through a Senior Research Fellowship.


Jean-François Coeurjolly
email: francois.coeurjolly@upmf-grenoble.fr

Pierre-Olivier Amblard

Sophie Achard

Wavelet analysis of the multivariate fractional Brownian motion

Keywords: multivariate fractional Brownian motion, wavelet analysis, cross-correlation, cross-spectrum

The work developed in the paper concerns the multivariate fractional Brownian motion (mfBm) viewed through the lens of the wavelet transform. After recalling some basic properties on the mfBm, we calculate the correlation structure of its wavelet transform. We particularly study the asymptotic behavior of the correlation, showing that if the analyzing wavelet has a sufficient number of null first order moments, the decomposition eliminates any possible long-range (inter)dependence. The cross-spectral density is also considered in a second part. Its existence is proved and its evaluation is performed using a von Bahr-Essen like representation of the function sign(t)|t| α . The behavior of the cross-spectral density of the wavelet field at the zero frequency is also developed and confirms the results provided by the asymptotic analysis of the correlation.

Motivations and overlook

The fractional Brownian (fBm) motion developed by [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF] has been extensively studied as the archetypal model of fractal signals. Many extensions have been proposed, trying to keep the simplicity of its definition while modeling more complex phenomena. For example, time-dependent Hurst exponent or d-dimensional extensions have been introduced which have respectively led to the multifractional Brownian motion [START_REF] Peltier | Multifractional Brownian motion: definition and preliminary results[END_REF] and the fractional Brownian sheet [START_REF] Ayache | Drap brownien fractionnaire[END_REF]. Another extension consists in defining multivariate fractal processes. This extension is needed by many applications ranging from economy to physics, passing by biology and neuroscience [START_REF] Gil-Alana | A fractional multivariate long memory model for the US and the Canadian real output[END_REF][START_REF] Arianos | Cross-correlation of long range correlated series[END_REF][START_REF] Achard | Fractal connectivity of long-memory networks[END_REF][START_REF] Achard | A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs[END_REF]. In all these disciplines, many modern sensing approaches allow to measure instantaneously different variables from complex phenomena. The need of multivariate signals models is crucial in order to model and understand these phenomena.

A multivariate extension of the fBm has been proposed recently in a very general setting by [START_REF] Didier | Integral representations of operator fractional Brownian motions[END_REF] with the help of operator self-similarity. The operator fractional Brownian motion is an operator self-similar Gaussian process with stationary increments. When the operator is diagonal, it is called the multivariate fractional Brownian motion (mfBm). We have particularly studied this diagonal case [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF], elaborating on the work of [START_REF] Lavancier | Covariance function of vector self-similar processes[END_REF]. In these works, the correlation structure of the mfBm has been studied. The increments process has also been studied, and we have exhibited its correlation and spectral structures, showing the possible existence of long-range dependence in correlation as well as in cross-correlation between components.

In this paper, we pursue the study by wavelet analyzing the multivariate process. It is now wellaccepted that wavelet analysis is the right framework to deal with monovariate fractal signals with stationary increments [START_REF] Flandrin | On the spectrum of fractional Brownian motions[END_REF][START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF][START_REF] Tewfik | Correlation structure of the discrete wavelet coefficients of fractional Brownian motion[END_REF][START_REF] Wornell | A Karhunen-Loève-like expansion for 1/f processes via wavelets[END_REF][START_REF] Faÿ | Estimators of long-memory: Fourier versus wavelets[END_REF][START_REF] Bardet | Wavelet estimator of long-range dependent processes[END_REF]. Wavelet transforms or decomposition provide a regularized differentiation of the processes, have a filter bank structure in perfect adequacy with 1/f type of spectral behavior, may eliminate long-range dependence properties if the analyzing wavelet is properly chosen. Thus, studying multivariate fractal signals through the lens of the wavelet zoom is indeed natural, and we expect that it will be useful as well in revealing the interaction structure between the components of the mfBm.

We thus concentrate on the correlation and on the spectral structure of the wavelet transform of the mfBm. The principal result of the paper is the explicit form of the second order statistics of this multivariate Gaussian random field. We study the asymptotic behavior of the cross-correlation function of the wavelet, and this allows us to prove that choosing a wavelet with at least two null first order moments eliminates any possible long dependence in the correlations. We prove the existence and calculate the cross-spectral density of the wavelet field. We thus extend the result of [START_REF] Kato | On the spectral density of the wavelet transform of fractional Brownian motions[END_REF] providing the existence of the spectral density of the wavelet transform of the fBm. The proof uses a generalization of the von Bahr-Essen representation of |t| α (used by [START_REF] Kato | On the spectral density of the wavelet transform of fractional Brownian motions[END_REF]) to the function sign(t)|t| α .

We also provide the behavior of the density at the zero frequency, corroborating the asymptotic result obtained in the time domain. We stress on the fact that our asymptotic results do not impose that the wavelet is a real function, nor that it has a compact support.

The paper is organized as follows. In order to have a self-contained exposition, we recall some basics definition and results on the mfBm in Section 1. In Section 2, we set the wavelet analysis and look at some self-similarity properties inherited from the process. The full correlation structure is developed in Section 3 whereas the spectral counterpart is in Section 4. The needed generalized von Bahr-Essen representation is proved in the last section.

Some facts on the multivariate fractional Brownian motion

The p dimensional multivariate fractional Brownian motion (mfBm) x(t) is defined as the only Gaussian process having stationary increments and having components jointly self-similar with parameters (H 1 , . . . , H p ) ∈ (0, 1) p . The self-similarity property can be stated as follows: for any real λ > 0,

x(λt) fidi = λ H x(t) where H = diag(H 1 , . . . , H p ) and λ H is intended in the matrix sense. The notation fidi = stands for equality of all the finite-dimensional probability distributions.

The cross-covariance structure induced by the multivariate self-similarity property and the stationarity of the increments has been first studied by [START_REF] Lavancier | Covariance function of vector self-similar processes[END_REF], Theorem 2.1, without having recourse to the Gaussian assumption. [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF] have parameterized this covariance structure in a more simple way as follows.

Proposition 1 (Proposition 3, [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF]) Let j, k ∈ {1, . . . , p}, j = k, then there exists σ j > 0, (ρ jk , η jk ) ∈ [-1, 1] × R satisfying ρ jk = ρ kj and η jk = -η kj , such that

r jk (s, t) := E[x j (s)x k (t)] = σ j σ k 2 {w jk (-s) + w jk (t) -w jk (t -s)} , (1)
where the function w jk (h) is defined by

w jk (h) =    (ρ jk -η jk sign(h))|h| Hj +H k if H j + H k = 1, ρ jk |h| + η jk h log |h| if H j + H k = 1. ( 2
)
This result is also valid in the case j = k when setting ρ jj = 1 and η jj = 0 in (2); we thus recover the covariance structure of a monovariate fBm. The parameter σ 2 j is the variance of a fBm at time 1, Var(x j (1)), whereas ρ jk represents the instantaneous correlation between components j and k at time 1, i.e. E[x j (1)x k (1)]. The antisymmetric parameter η jk is related to the time-reversibility property of the multivariate process. Indeed, the mfBm is time reversible, i.e. x(t) = x(-t) in distribution for every t, if and only if η jk = 0 for all j, k [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF].

To ensure that the matrix given by ( 1) is the cross-covariance matrix of a process, the constraints imposed on ρ jk and η jk are not sufficient. A necessary and sufficient condition, proved by [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF], corresponds to the positive-definiteness of the Hermitian matrix with entries Γ(H j + H k + 1) × ξ jk where ξ jk is defined by

ξ jk =    ρ jk sin π 2 (H j + H k ) -iη jk cos π 2 (H j + H k ) if H j + H k = 1 ρ jk -i π 2 η jk if H j + H k = 1, (3)
where i = √ -1. For example, when p = 2 and η jk = 0 no condition is required for the correlation ρ 12 but when H 1 = H 2 , and when H 1 = 0.1 and H 2 = 0.2 the correlation ρ 12 cannot exceed 0.514, see [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF] for more discussion. The problem of simulation of such a process has been investigated in [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF] using the algorithm developed by [START_REF] Chan | Simulation of stationary Gaussian vector fields[END_REF]. The covariance structure of the increment process (at lag 1) can be easily deduced from (1). When j = k, we obviously recover the covariance of the fractional Gaussian noise and the classical property that this process is short-memory when H j ≤ 1/2 and long-memory when H j > 1/2. In the multivariate case, long-range (interdependence) may also appear in the cross-covariance. Indeed [START_REF] Amblard | Basic properties of the multivariate fractional Brownian motion[END_REF] proved that for all j = k, the cross-covariance behaves asymptotically as |h| Hj +H k -2 (up to a constant) meaning that the long-memory property arises as soon as H j +H k ≥ 1 which can appear in three different situations:

H j = H k = 1/2, H j < 1/2 and H k > 1 -H j or H j > 1/2 and H k > 1/2.
In those cases, some troubles may appear when it comes to infer parameters of the models from data. Indeed, longrange dependence may lead to very slow convergence of estimators. As already observed in many works, e.g. [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF][START_REF] Veitch | Wavelet-based joint estimate of the long-range dependence parameters[END_REF][START_REF] Bardet | Wavelet estimator of long-range dependent processes[END_REF][START_REF] Coeurjolly | Estimating the parameters of a fractional Brownian motion by discrete variations of its sample paths[END_REF], having recourse to wavelet types of transformation is an elegant way to overcome the problem. Indeed, using wavelet types of transformation with a correctly chosen filter allows to extract the stationary part from the fBm and allows to "whiten" the increments. We describe such an approach in the following section.

Wavelet Analysis: definition, stationarity and self-similarity

The use of wavelet analysis in the understanding of the monovariate fractional Brownian motion, and more generally for the study of fractal signals, goes back to the early works of [START_REF] Flandrin | On the spectrum of fractional Brownian motions[END_REF][START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF], [START_REF] Tewfik | Correlation structure of the discrete wavelet coefficients of fractional Brownian motion[END_REF], [START_REF] Wornell | A Karhunen-Loève-like expansion for 1/f processes via wavelets[END_REF] to cite some but a few.

The aim is now to analyze the multivariate fractional Brownian motion through the lens of the wavelet transform. We use the continuous wavelet transform here, but a similar analysis could be performed in the multiresolution framework using orthonormal wavelet bases. We will consider complex valued wavelets, not necessarily in the Hardy class, not necessarily with compact support. The hypothesis we impose on the wavelets will be detailed when needed.

Definition and stationarity

Let ψ be a complex wavelet function, let a > 0 and b ∈ R and consider ψ ab (.

) = a -1/2 ψ((. -b)/a). Let d j a,b := x j ψ a,b L 2 = R x j (t)ψ a,b (t)dt (4)
the wavelet transform of the jth component of a multivariate fractional Brownian motion. ψ denotes the complex conjugate of ψ. In this section, we assume that conditions [C1] and [C2(2)] are satisfied, where:

[C1] Admissibility condition:

ψ(t) ∈ L 2 and | ψ(ω)| 2 /|ω| ∈ L 1
, where ψ is the Fourier transform of ψ.

[C2(K)] t m ψ(t) ∈ L 1 for m = 0, 1, . . . , K.

Condition [C1] ensures that ψ(0) = 0 and that R ψ(t)dt = 0. We note, as [START_REF] Kato | On the spectral density of the wavelet transform of fractional Brownian motions[END_REF], that under condition [C2(1)], the integral ( 4) is well-defined as a sample path integral and is a second-order random variable. This follows, since under [C2(1)] we have R |s| H |ψ a,b (s)|ds < +∞, ∀H ∈ (0, 1).

The aim of this section is to focus on the correlation between the wavelet transforms (at different scales and different times) of two components j and k of the multivariate fractional Brownian motion.

The wavelet transform is a random field. It is clearly zero mean and Gaussian. We have for a 1 , a 2 > 0 and b, h ∈ R

E[d j a1,b+h d k a2,b ] = R 2 r jk (t 1 , t 2 )ψ a1,b+h (t 1 )ψ a2,b (t 2 )dt 1 dt 2 .
Under [C1], and from (2) the last expression reduces to

E[d j a1,b+h d k a2,b ] = - σ j σ k 2 R 2 w jk (t 2 -t 1 )ψ a1,b+h (t 1 )ψ a2,b (t 2 )dt 1 dt 2 . Let Γ ψ (v) := R ψ a1,b+h (u)ψ a2,b (u + v
)du be the correlation function between the two wavelets ψ a1,b+h and ψ a2,b . Then we have

E[d j a1,b+h d k a2,b ] = - σ j σ k 2 R w jk (v)Γ ψ (v)dv. (5)
Note that [C2(2)] implies that for all the values of H j and

H k , R |w jk (v)| |Γ ψ (v)| dv < +∞.
With two changes of variables, this may also be rewritten as

E[d j a1,b+h d k a2,b ] = - σ j σ k 2 √ a 1 a 2 × R 2 w jk (a 2 t 2 -a 1 t 1 -h)ψ(t 1 )ψ(t 2 )dt 1 dt 2 . ( 6
)
If we interpret for fixed parameters a 1 and a 2 , the quantity E[d j a2,b+h d k a2,b ] as the cross-correlation between two signals, we observe that it depends only on the difference between the times at which it is evaluated (i.e. h). With the fact that the wavelet transform is a zero mean and Gaussian field, we conclude that d j a1,. and d k a2,. are jointly stationary signals.

Self-similarity type property of the cross-wavelet transform

The variance of the wavelet transform at similar scales for the fractional Brownian motion with Hurst parameter H exhibits some self-similarity. Indeed, it is proved in [START_REF] Flandrin | On the spectrum of fractional Brownian motions[END_REF] for example that for

all b V ar(d j a,b ) = a 2H+1 × - σ 2 2 R 2 |t 2 -t 1 | 2H ψ(t 1 )ψ(t 2 )dt 1 dt 2 .
We note here that the same behavior holds for the cross-wavelet variance.

Proposition 1 Under the assumptions [C1] and [C2(2)]

, let h = 0 and fix a 1 = a 2 = a > 0. Then,

E[d j a,b d k a,b ] = a Hj +H k +1 - σ j σ k 2 z jk and Corr[d j a,b , d k a,b ] = z jk √ z jj z kk , (7)
where

z jk := R 2 w jk (t 2 -t 1 )ψ(t 1 )ψ(t 2 )dt 1 dt 2 .
Proof.

Consider Equation ( 6). The result is obvious when H j + H k = 1 since for any a > 0,

w jk (av) = a Hj +H k w jk (v). Now, when H j + H k = 1, the result comes from Condition [C1] ensuring that R 2 η jk × a(t 2 -t 1 ) log(a)ψ(t 1 )ψ(t 2 )dt 1 dt 2 = 0.
Let us observe that the instantaneous cross-wavelet correlation is independent of the scale.

4 Cross-correlation structure of the wavelet transform of the mfBm For fixed scales, a 1 , a 2 , we now specify the behavior of the cross-wavelet covariance (or correlation) as |h| → +∞. In particular, our aim is to exhibit the influence of the number of vanishing moments of the wavelet function on the asymptotic cross-wavelet covariance. Such a result needs the following assumption:

[C3] The wavelet function has M ≥ 1 vanishing moments that is We may now derive our result obtained as |h| → +∞. Let us first recall Landau notation: for two functions f (h) and g(h) defined on R, we denote by f (h)∼g(h) as |h| → +∞ (resp. f (h)=o(g(h)) and )] and [C3] hold, then as |h| → +∞, we have

f (h)=O(g(h))) if lim |h|→+∞ f (h)/g(h) = 1 (resp. lim |h|→+∞ f (h)/g(h) = 0 and f (h)/g(h) is bounded for all h). Theorem 2 Assume [C1], [C2(2M+1
E[d j a1,b+h d k a2,b ] ∼ - σ j σ k 2 κ(ψ, M ) τ jk (h) |h| Hj +H k -2M
where κ(ψ, M ) :

= 2M M (a 1 a 2 ) M t M ψ(t)dt 2 and τ jk (h) =    (ρ jk + η jk sign(h)) Hj +H k 2M if H j + H k = 1 - η jk ×sign(h) 2M(2M-1) if H j + H k = 1. ( 8
)
We notice that the equivalence stated has a meaning as soon as ρ jk + η jk sign(h) = 0 when H j + H k = 1 and as soon as η jk = 0 when H j + H k = 1. In the opposite cases, a careful look at the proof shows that the equivalence can be replaced by an upper-bound or more precisely

E[d j a1,b+h d k a2,b ] = o(|h| 1-2M ).
Similarly to the fractional Brownian motion, Theorem 2 asserts that the higher M , the less correlated the wavelet transforms of the components j and k of the multivariate fractional Brownian motion.

This has many implications. In particular, this suggests that estimating the instantaneous cross-wavelet correlation at a scale a may be efficiently done by using the empirical correlation since at scale a, d j a,b+h and d k a,b are not too much correlated if M is large.

Proof. The proof is split into two cases. Before this, we denote by

D h := {(t 1 , t 2 ) ∈ R 2 : |a 2 t 2 -a 1 t 1 | < |h| 2
} and we note in particular that

∀(t 1 , t 2 ) ∈ D h a2t2-a1t1 h ≤ 1 2 < 1 and sign(a 2 t 2 -a 1 t 1 -h) = -sign(h), ∀(t 1 , t 2 ) ∈ R 2 \ D h h a2t2-a1t1 ≤ 2. Case 1. α := H j + H k = 1.
We assume here that

ρ jk + η jk sign(h) = 0. Let us write E[d j a,b+h d k a,b ] = - σj σ k 2 √ a 1 a 2 × T with T := R 2 (ρ jk -η jk sign(a 2 t 2 -a 1 t 1 -h)|a 2 t 2 -a 1 t 1 -h| α ψ(t 1 )ψ(t 2 )dt 1 dt 2 = T 1 + T 2 ,
and where T 1 (resp. T 2 ) corresponds to the integral on D h (resp. R 2 \ D h ). Let us first prove that

|h| 2M-α T 2 → 0 as |h| → +∞. Denoting c ∨ = |ρ jk | + |η jk |, we have (since 2M -α > 0) |h| 2M-α |T 2 | ≤ c ∨ R 2 \D h |a 2 t 2 -a 1 t 1 | α |h| 2M-α 1 + h a 2 t 2 -a 1 t 1 α |ψ(t 1 )||ψ(t 2 )|dt 1 dt 2 ≤ 2 2M-α 3 α c ∨ R 2 \D h (a 2 t 2 -a 1 t 1 ) 2M |ψ(t 1 )||ψ(t 2 )|dt 1 dt 2 .
The result is then obtained by using assumption [C2(2M)] and the dominated convergence theorem. Now, within the domain D h , one may use the series expansion of (1 + x) α (for |x| < 1).

T 1 = |h| α D h (ρ jk -η jk sign(a 2 t 2 -a 1 t 1 -h)) 1 - a 2 t 2 -a 1 t 1 h α ψ(t 1 )ψ(t 2 )dt 1 dt 2 = |h| α (ρ jk + η jk sign(h)) D h   ℓ≥0 (-1) ℓ α ℓ a 2 t 2 -a 1 t 1 h ℓ   ψ(t 1 )ψ(t 2 )dt 1 dt 2 ,
where α ℓ denotes the binomial coefficient (α)(α -1) . . . (α -ℓ + 1)/ℓ!. Decompose T 1 into three terms (denoted by T ′ 1 , T ′ 2 and T ′ 3 ) corresponding to the 2M first terms of the series, the (2M + 1)th term (ℓ = 2M ) and the remainder terms. Then,

T ′ 1 = |h| α (ρ jk + η jk sign(h)) 2M-1 ℓ=0 (-1) ℓ h -ℓ α ℓ D h (a 2 t 2 -a 1 t 1 ) ℓ ψ(t 1 )ψ(t 2 )dt 1 dt 2 .
Under Assumption [C3], ψ has M vanishing moments and therefore the previous expression reduces to

T ′ 1 = -|h| α (ρ jk + η jk sign(h)) 2M-1 ℓ=0 (-1) ℓ h -ℓ α ℓ R 2 \D h (a 2 t 2 -a 1 t 1 ) ℓ ψ(t 1 )ψ(t 2 )dt 1 dt 2 . Now, |h| 2M-α |T ′ 1 | ≤ c ∨ 2M-1 ℓ=0 α ℓ R 2 \D h 2 2M-ℓ (a 2 t 2 -a 1 t 1 ) 2M |ψ(t 1 )||ψ(t 2 )|dt 1 dt 2 .
Assumption [C2(2M)] and the dominated convergence theorem may be combined to prove that |h| 2M-α T ′ 1 → 0 as |h| → +∞. The term T ′ 2 is defined as

T ′ 2 := |h| α-2M (ρ jk + η jk sign(h)) α 2M D h (a 2 t 2 -a 1 t 1 ) 2M ψ(t 1 )ψ(t 2 )dt 1 dt 2 .
As previously we obtain

|h| 2M-α T ′ 2 ρ jk + η jk sign(h) → α 2M R 2 (a 2 t 2 -a 1 t 1 ) 2M ψ(t 1 )ψ(t 2 )dt 1 dt 2 = α 2M 2M M (a 1 a 2 ) M t M ψ(t)dt 2 = α 2M κ(ψ, M ). Since T = T 1 + T 2 = T ′ 1 + T ′ 2 + T ′ 3 + T 2 the proof will be completed if we manage to prove that |h| 2M-α T ′ 3 → 0. Let us write |h| 2M-α T ′ 3 = h 2M (ρ jk + η jk sign(h)) D h ℓ≥2M+1 (-1) ℓ α ℓ a 2 t 2 -a 1 t 1 h ℓ ψ(t 1 )ψ(t 2 )dt 1 dt 2 = ρ jk + η jk sign(h) h × D h (a 2 t 2 -a 1 t 1 ) 2M+1   ℓ≥0 (-1) ℓ+1 α ℓ + 2M + 1 a 2 t 2 -a 1 t 1 h ℓ   ψ(t 1 )ψ(t 2 )dt 1 dt 2 .
The binomial coefficient appearing in the last equation satisfies, with

ℓ ′ = ℓ + 2M + 1 α ℓ ′ = α(α -1) • • • (α -ℓ ′ + 1) ℓ ′ ! ≤ 2(2 -α) • • • (ℓ ′ -1 -α) ℓ ′ ! since α ≤ 2 ≤ 2(ℓ ′ -1)! ℓ ′ ! = 2 ℓ ′ ≤ 2 ℓ .
Recall that in D h we have a 2 t 2 -a 1 t 1 /|h| ≤ 1/2. The series in the previous integral then satisfies

ℓ≥0 α ℓ + 2M + 1 a 2 t 2 -a 1 t 1 h ℓ ≤ α 2M + 1 + ℓ≥1 α ℓ + 2M + 1 a 2 t 2 -a 1 t 1 h ℓ ≤ 2 2M + 1 + ℓ≥1 2 ℓ 2 -ℓ = 2 2M + 1 + 2 log(2) =: C M .
Thus we obtain

|h| 2M-α T ′ 3 ≤ C M c ∨ |h| R 2 a 2 t 2 -a 1 t 1 2M+1 |ψ(t 1 )||ψ(t 2 )|dt 1 dt 2 . Since by Assumption [C2(2M+1)], t 2M+1 ψ(t) ∈ L 1 , we have |h| 2M-α T ′ 3 = O(|h| -1
), whence the result.

Case 2.

H j + H k = 1.
We assume here that η jk = 0. We take the same notation as previously. We first note that, under

[C1], the term T can be rewritten as

T = R 2 ρ jk |a 2 t 2 -a 1 t 1 -h| + η jk (a 2 t 2 -a 1 t 1 -h) log 1 - a 2 t 2 -a 1 t 1 h ψ(t 1 )ψ(t 2 )dt 1 dt 2 .
We decompose T in T 1 + T 2 (as done in case 1). The proof that |h| 2M-1 T 2 → 0 as |h| → +∞ follows similar arguments as in the case 1 and is therefore omitted. Now, the term T 1 can be rewritten as

T 1 = ρ jk |h| D h 1 - a 2 t 2 -a 1 t 1 h ψ(t 1 )ψ(t 2 )dt 1 dt 2 -η jk h D h 1 - a 2 t 2 -a 1 t 1 h log 1 - a 2 t 2 -a 1 t 1 h ψ(t 1 )ψ(t 2 )dt 1 dt 2 .
Denote by T1 and T2 these two terms. Assumption [C1] leads to

T1 = -ρ jk |h| R 2 \D h 1 - a 2 t 2 -a 1 t 1 h ψ(t 1 )ψ(t 2 )dt 1 dt 2 = -ρ jk R 2 \D h |a 2 t 2 -a 1 t 1 | 1 - h a 2 t 2 -a 1 t 1 dt 1 dt 2 .
Then, we assert that

|h| 2M-1 | T1 | ≤ 2 2M-1 3|ρ jk | R 2 \D h a 2 t 2 -a 1 t 1 2M |ψ(t 1 )||ψ(t 2 )|dt 1 dt 2 → 0
as |h| → +∞. For the term T2 , we may use the series expansion of log(1 + x) (for |x| < 1). We omit the details and leave the reader to verify that as |h| → +∞

T2 ∼ η jk h R 2 1 - a 2 t 2 -a 1 t 1 h -1 2M -1 a 2 t 2 -a 1 t 1 h 2M-1 -1 2M a 2 t 2 -a 1 t 1 h 2M ψ(t 1 )ψ(t 2 )dt 1 dt 2 ∼ -h 1-2M × η jk 2M (2M -1) 2M M (a 1 a 2 ) M t M ψ(t)dt 2 . Hence, T ∼ |h| 1-2M × - η jk ×sign(h) 2M(2M-1) κ(ψ, M ).
In this proof, Fubini's theorem and interchanges of integrals and (in)finite sums are widely used. All of these are justified by the absolute convergence of the different series related to the expansions of (1 + x) α

or log(1 + x) for |x| < 1 and Assumption [C2(2M+1)].

5 Cross-spectral density of the wavelet transform of the mfBm

In the case of the fBm, the expression of the spectral density of the wavelet transform was provided by [START_REF] Flandrin | On the spectrum of fractional Brownian motions[END_REF][START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF]. A rigorous proof of the existence of this spectral density in the L 1 sense was obtained by [START_REF] Kato | On the spectral density of the wavelet transform of fractional Brownian motions[END_REF]. On the basis of this work, our ambition is to provide the cross-spectral density between wavelet transforms (at different scales) of components j and k of the multivariate fractional Brownian motion. The idea is to obtain the following spectral representation for the cross-correlation

E[d j a1,b+h d k a2,b ] = 1 2π R S jk a1,a2 (ω)
e iωt dω.

Theorem 3 Under Assumptions [C1], [C2(M)] and [C3]

(with M ≥ 2), we derive the following assertions.

(i) The cross-spectral density of the wavelet transforms of two components j and k exists and is given by

S jk a1,a2 (ω) = √ a 1 a 2 σ j σ k Γ(H j + H k + 1) ζ jk (ω) ψ(a 1 ω) ψ(a 2 ω) |ω| Hj +H k +1 (9)
where

ζ jk (ω) =    ρ jk sin π 2 (H j + H k ) + i η jk cos π 2 (H j + H k ) sign(ω) if H j + H k = 1 ρ jk + i π 2 η jk sign(ω) if H j + H k = 1.
(ii) We have for both cases, as ω → 0

S jk a1,a2 (ω) ∼ (a 1 a 2 ) M+1/2 σ j σ k Γ(H j + H k + 1) | ψ (M) (0)| 2 |ζ jk (ω)| |ω| 2M-1-α .
(iii) Moreover, the coherence function between the two components j and k satisfies:

C jk a1,a2 (ω) := S jk a1,a2 (ω) 2 S jj a1,a1 (ω)S kk a2,a2 (ω) = |ζ jk (ω)| 2 Γ(H j + H k + 1) 2 Γ(2H j + 1)Γ(2H k + 1) × ψ(a 1 ω) ψ(a 2 ω) ψ(a 1 ω) ψ(a 2 ω) . ( 10
)
Before writing down the proof, let us give some comments.

1. Item (ii) in Theorem 3 is the spectral analogue of Theorem 2. Indeed the behavior ot the crosscorrelation at infinite lags is linked to the behavior of its Fourier transform at the zero frequency.

We recover the fact that as soon as M > H j + H k + 1/2, the long-range interdependence is destroyed. The divergence of |ω| -1-Hj +H k is compensated by the rapid decrease to zero of the Fourier transform of the wavelet.

2. The interpretation of the coherence ( 10) is difficult here. Indeed, it is complex valued, a property which is not natural for a coherence. This comes from the fact that the quantities S jj a1,a2 (ω) are not power spectral densities but cross-spectral densities (cross-spectral density between two different scales of the wavelet transform of one signal). Thus, to interpret correctly the coherence, we should look at one scale only, in which case we recover the coherence evaluated in the usual spectral domain.

And this result is logical since the usual coherence is independent of the frequency.

3. Setting a 1 = a 2 and j = k in the expression of the cross-spectral density, we recover the usual result of the power spectral density at one scale of the wavelet transform of a scalar fBm. The proof proposed here is a natural extension of the proof found in [START_REF] Kato | On the spectral density of the wavelet transform of fractional Brownian motions[END_REF].

4. The derivation of the analytic form of the cross-spectral density is easy if we use generalized functions (or Schwartz distributions). Indeed, from [START_REF] Gel | Properties and Operations[END_REF], we know that the Fourier transforms of |v| α and |v| α sign(v) are respectively given by -2Γ(α + 1)|ω| -α-1 sin(πα/2) and -2iΓ(α + 1)|ω| -α-1 sign(ω) sin(πα/2). Then, in the proof below, the calculation of T = R w jk (v)Γ ψ (v)dv can be done using Parseval equality. However, the theoretical background required and hidden in the calculation is far more involved than the basics we have used in the proof (see [START_REF] Gel | Properties and Operations[END_REF]).

Proof. (i) We recall that under [C1] and [C2(2)], Equation ( 5) holds, that is E[d j a1,b1 d k a2,b2 ] = -σj σ k 2 T with T := R w jk (v)Γ ψ (v)dv. Furthermore, note that the Fourier transforms of ψ a,b and Γ ψ (v) exist and are equal respectively to √ a ψ(aω)e -iωb and to

q(ω) := Γ ψ (ω) = R Γ ψ (v)e -iωv dv = √ a 1 a 2 ψ(a 1 ω) ψ(a 2 ω)e iωh . (11)
Now, let us split the proof into two cases.

Case 1. α := H j + H k = 1.

When j = k, at this step, [START_REF] Kato | On the spectral density of the wavelet transform of fractional Brownian motions[END_REF] have used the representation of |v| α obtained by von [START_REF] Bahr | Inequalities for the r-th absolute moment of a sum of random variables, 1 ≤ r ≤ 2[END_REF]. We have obtained a similar representation for the function sign(v)|v| α for α ∈ (0, 2) \ {1} (see Equations ( 13) and ( 14) in Lemma 4). We have by Fubini's theorem and under

Assumption [C3] (with M ≥ 2). T = R (ρ jk -η jk sign(v))|v| α Γ ψ (v)dv = Γ(α + 1) π R |ω| -α-1 R ρ jk sin(πα/2)(1 -cos(ωv)) -η jk cos(πα/2)sign(ω)(sin(ωv) -g α (ωv)) × Γ ψ (v) dv dω = Γ(α + 1) π R |ω| -α-1 -ρ jk sin(πα/2) q(-ω) + q(ω) 2 -η jk cos(πα/2)sign(ω) q(-ω) -q(ω) 2i dω = - Γ(α + 1) π R |ω| -α-1 ρ jk sin(πα/2) + i η jk cos(πα/2)sign(ω) q(ω)dω = - Γ(α + 1) π R |ω| -α-1 ζ jk (ω)q(ω)dω.
Note that the condition M ≥ 2 is required for α > 1. For α < 1, M ≥ 1 is a sufficient condition. These conditions allow us to show that the contributions Γ ψ (v)dv and R g α (ωv)Γ ψ (v)dv are equal to zero. Now, using (11) we obtain

E[d j a1,b1 d k a2,b2 ] = √ a 1 a 2 σ j σ k Γ(α + 1) × 1 2π R |ω| -α-1 ζ jk (ω) ψ(a 1 ω) ψ(a 2 ω) =:P (ω)
e iωh dω.

By using Bochner's theorem, the proof will be done, if one proves that the function P (•) is integrable.

Let us prove this last assertion. Under [C2(M)], t k ψ(t) ∈ L 1 for k = 0, . . . , M . Therefore, ψ is a M times continuous and differentiable function. Using a Taylor expansion

ψ(ω) = M-1 k=0 ω k ψ (k) (ω) + ω M ψ (M) ( ω) = ω M ψ (M) ( ω), with ω ∈ [0 ∧ ω, 0 ∨ ω],
under [C2(M)]. And since ψ (M) is continuous at zero, ψ(ω) ∼ ω M ψ (M) (0) as ω → 0. Then as ω → 0:

P (ω) ∼ ζ jk (ω)|ω| 2M-1-α (a 1 a 2 ) M | ψ (M) (0)| 2 . ( 12
)
As a consequence, for M ≥ 2, P is continuous at zero and lim ω→0 ± P (ω) = 0. Therefore for ε > 0, P is integrable on the interval [-ε, ε] as a continuous function on this interval. Finally (with

c ∨ := |ρ jk |+|η kj |), |ω|≥ε |P (ω)| ≤ c ∨ a α 1 |ω|≥a1ε | ψ(ω)| 2 |ω| α+1 dω 1/2 a α 2 |ω|≥a2ε | ψ(ω)| 2 |ω| α+1 dω 1/2 ≤ c ∨ ε α R | ψ(ω)| 2 |ω| dω < +∞, under [C1]
. Hence, P (•) ∈ L 1 and Bochner's theorem may be applied.

Case 2. H j + H k = 1.

We start with the representation of v log |v| given by ( 16).

w jk (v) = ρ jk |v| + η jk v log |v| = lim α→1 -ρ jk |v| α + η jk v log |v| = lim α→1 - 1 2π R 2ρ jk (1 -cos(ωv)) -πη jk sign(ω) sin(ωv) |ω| α+1 dω.
Now, we derive the computation of the term T := R w jk (v)Γ ψ (-v)dv, similarly as the previous case.

Using dominated convergence theorem and Fubini's theorem,

T = 1 2π R lim α→1 -R 2ρ jk (1 -cos(ωv)) -πη jk sign(ω) sin(ωv) |ω| α+1 dω Γ ψ (v)dv = 1 2π lim α→1 -R R 2ρ jk (1 -cos(ωv)) -πη jk sign(ω) sin(ωv) |ω| α+1 Γ ψ (v)dv dω = 1 2π lim α→1 -R -2ρ jk q(-ω) + q(ω) 2 -πη jk sign(ω) q(-ω) -q(ω) 2i |ω| -α-1 dω = - 1 2π lim α→1 -R 2ρ jk + iπη jk sign(ω) |ω| α+1 q(ω)dω = - 1 2π lim α→1 -R
|ω| -α-1 (2ρ jk + iπη jk sign(ω)) ψ(a 1 ω) ψ(a 2 ω)e iωh dω.

From ( 12), |ω| -α-1 ψ(a 1 ω) ψ(a 2 ω) is an integrable function for all α ∈ (0, 2). Therefore, the integral and the limit may be interchanged. Therefore, we obtain The following lemma provides a similar representation for sign(v)|v| α , v α dω.

+ = v α 1 R + (v) and v α -= (-v) α 1 R -(v).
The representations of v α + and v α -are obtained from ( 13) and ( 14) noticing that

FigureFigure 1 :

 1 Figure 1: Examples of discretized sample paths of a time reversible (η jk = 0) mfBm of length n = 1024, with p = 20 components. The Hurst exponents are equally spaced in [0.3, 0.4] (left plot), [0.6, 0.7] (middle plot) and [0.4, 0.8] (right plot). The correlation parameters are set to 0.7 (left and middle plot) and to 0.3 (right plot). The components are shifted artificially in the left plot for the sake of visibility.

R

  t m ψ(t)dt = 0 for m = 0, . . . , M -1 and R t M ψ(t)dt = 0.

ρ

  jk + i π 2 η jk sign(ω) |ω| 2 ψ(a 1 ω) ψ(a 2 ω)e iωh dω, and Bochner's theorem can be applied.(ii) is derived from (12). 6 Bahr and Essen type representations for the functions sign(v)|v| α , v α + and v α -In 1965, von Bahr and Essen have obtained the following representation theorem for |v| α for α ∈ (0, 2):

Lemma 4

 4 Let α ∈ (0, 2) \ {1} and let g α : R → R the function which equals zero when α ∈ (0, 1) and which is the identity function when α ∈ (1, 2), then we have sign(v)|v| α = Γ(α + 1) cos(πα/2) π R sign(ω) (sin(ωv) -g α (ωv)) cos(ωv)) + cos π α 2 sign(ω) (sin(ωv) -g α (ωv)) cos(ωv)) -cos π α 2 sign(ω) (sin(ωv) -g(ωv)) |ω| α+1

  α -sign(v)|v| α ) .Proof. Let α ∈ (0, 1), then from (13) and properties of the function Γ1 α + 1 |v| α+1 = Γ(α + 1) π cos(πα/2) R 1 -cos(ωv) |ω| α+2 dω.Since R |ω| -α-1 | sin(ωv)| < +∞ for α ∈ (0, 1), we can differentiate this integral with respect to vSince R |ω| -α-1 | sin(ωv) -ωv|dω < +∞ for α ∈ (1, 2), we can take the primitive of the last equation )(sin(ωv) -ωv) |ω| α+1 dω, which ends the proof.Let α ∈ (0, 1), then by differentiating (15) with respect to α and taking the limit as α → 1 -, we may obtain sign(h)|h| log |h| = h log |h| = lim α→1 -