Charles-Gérard Lucas

Gustavo Didier
email: gdidier@tulane.edu

Herwig Wendt
email: herwig.wendt@irit.fr.

Ieee Senior

Patrice Abry
email: patrice.abry@ens-lyon.fr

Multivariate selfsimilarity: Multiscale eigenstructures for selfsimilarity parameter estimation

Keywords: Multivariate selfsimilarity, selfsimilarity parameter vector, wavelet transform, multiscale eigen-structure

Scale-free dynamics, formalized by selfsimilarity, provides a versatile paradigm massively and ubiquitously used to model temporal dynamics in real-world data. However, its practical use has mostly remained univariate so far. By contrast, modern applications often demand multivariate data analysis. Accordingly, models for multivariate selfsimilarity were recently proposed. Nevertheless, they have remained rarely used in practice because of a lack of available robust estimation procedures for the vector of selfsimilarity parameters. Building upon recent mathematical developments, the present work puts forth an efficient estimation procedure based on the theoretical study of the multiscale eigenstructure of the wavelet spectrum of multivariate selfsimilar processes. The estimation performance is studied theoretically in the asymptotic limits of large scale and sample sizes, and computationally for finite-size samples. As a practical outcome, a fully operational and documented multivariate signal processing estimation toolbox is made freely available and is ready for practical use on real-world data. Its potential benefits are illustrated in epileptic seizure prediction from multi-channel EEG data.

I. INTRODUCTION

Context. The concept of scale-free dynamics results from a change in paradigm: instead of being driven by a few characteristic time scales, thus targeted in estimation, temporal dynamics results from the interactions of a large continuum of time scales, hence all of equal relevance. Scale-free dynamics has been massively used in the past decades to successfully model temporal activities in a large variety of data stemming from several real-world applications very different in nature, ranging from natural phenomena -physics (hydrodynamic turbulence [START_REF] Frisch | Turbulence: the Legacy of A.N. Kolmogorov[END_REF], astrophysics [START_REF] Davis | Cosmic microwave background probes models of inflation[END_REF]), geophysics (rainfalls, [START_REF] Hubert | Multifractal taming of extreme hydrometeorological events[END_REF], seismology [START_REF] Turcotte | Implications of chaos, scale-invariance, and fractal statistics in geology[END_REF]), biology (plant growths [START_REF] Palmer | Fractal geometry: a tool for describing spatial patterns of plant communities[END_REF], body rhythms [START_REF] Lopes | Fractal and multifractal analysis: A review[END_REF], heart rate [START_REF] Ivanov | Scale-invariant aspects of cardiac dynamics across sleep stages and circadian phases[END_REF], [START_REF] Nakamura | Multiscale analysis of intensive longitudinal biomedical signals and its clinical applications[END_REF], neurosciences [START_REF] He | Scale-free brain activity: past, present, and future[END_REF], [START_REF] La Rocca | Self-similarity and multifractality in human brain activity: a waveletbased analysis of scale-free brain dynamics[END_REF], genomics [START_REF] Takahashi | A fractal model of chromosomes and chromosomal dna replication[END_REF])to human activity -Internet traffic [START_REF] Fontugne | Scaling in internet traffic: a 14 year and 3 day longitudinal study, with multiscale analyses and random projections[END_REF], finance [START_REF] Mandelbrot | A multifractal walk down wall street[END_REF], city growth [START_REF] Lengyel | Roughness and intermittency within metropolitan regions-Application in three French conurbations[END_REF], art investigations [START_REF] Abry | Multiscale Anisotropic Texture Analysis and Classification of Photographic Prints: Art scholarship meets image processing algorithms[END_REF].

However, scale-free dynamics has so far remained restricted to univariate settings only, mostly due to a lack of multivariate models and estimation tools. By contrast, in most applications, one same system is monitored by several sensors, which Supported by PhD Grant DGA/AID (no 01D20019023) and ANR-18-CE45-0007 MUTATION and the Simons Foundation collaboration grant #714014.

naturally calls for the joint (or multivariate) analysis of a collection of time series. Multivariate selfsimilarity with the practical and accurate estimation of the vector of selfsimilarity parameters thus constitutes the heart of the present work. Related works. Fractional Brownian motion (fBm) [START_REF] Samorodnitsky | Stable non-Gaussian Random Processes: Stochastic Models with Infinite Variance[END_REF], [START_REF] Pipiras | Long-Range and Self-Similarity[END_REF], the only Gaussian, selfsimilar, stationary-increment process, has quasi-exclusively been used in practice as the reference model for scale-free dynamics. The selfsimilarity parameter H is the quantity of interest in applications. Its accurate estimation has thus received considerable attention (cf. [START_REF] Pipiras | Long-Range and Self-Similarity[END_REF], [18] for reviews). Notably, multiscale (wavelet) representations have proven to yield accurate, robust and efficient estimation procedures for H [START_REF] Wendt | [END_REF]: They rely on the crucial facts that the statistics of the wavelet coefficients behave as power laws with respect to scales, and that the exponents of these power laws are controlled by H. To account for the multivariate nature of data in applications, multivariate extensions of fBm were recently proposed based on gathering correlated fBm [20]. Multivariate multiscale-based representations were accordingly devised for the estimation of the resulting vector of selfsimilarity parameters. Nonetheless, these constructs were essentially based on component-wise univariate power law behavior [START_REF] Wendt | [END_REF], [START_REF] Coeurjolly | Wavelet the multivariate Brownian motion[END_REF]. More recently, a richer model for multivariate selfsimilarity, operator fractional Brownian motion (ofBm), was proposed. Its versatility stems from a major change of perspective: there is no longer a one-toone association either between each selfsimilarity parameter and each component, or between power laws and scale-free dynamics. Instead, the statistics of each component depend on the entire vector of selfsimilarity parameters [START_REF] Mason | Sample path properties of similar Gaussian fields[END_REF]- [START_REF]Exponents, symmetry groups and classification of operator frac-Brownian motions[END_REF]. Very recently, the principles of the estimation for the full vector of selfsimilarity parameters were theoretically discussed, based on the eigenvalues of multivariate wavelet spectra estimated at each scale [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF]Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]. This paved the way for preliminary explorations [START_REF] Lucas | Bootstrap for testing the equality of selfsimilarity exponents across multivariate time series[END_REF] that lead to the core contributions of the present work: The construction of practical statistical signal processing tools for multivariate selfsimilarity analysis. Goals, contributions and outline. The present work aims to define theoretically well-grounded and practically robust, efficient and operational estimation procedures for the vector of selfsimilarity parameters in a multivariate setting and to study theoretically and practically their performance. To that end, Section II recalls state-of-the-art multivariate selfsimilarity modeling and analysis. As a first contribution, a practically tractable yet richer model for multivariate selfsimilarity is devised in Section III, and then compared to ofBm; its multivariate wavelet analysis is detailed. As a second and major contribution, an estimation procedure for the vector of selfsimilarity parameters is devised, by exploiting a bias-variance tradeoff, and its asymptotic performance is theoretically studied (cf. Section IV). Finite-size estimation performance is studied empirically (from Monte Carlo simulations), together with the derivation of approximations to the covariance structure of the estimates, which is of critical importance for real-world applications (cf. Section V). The estimation performance of the proposed multivariate estimator is compared, in terms of biases, variances, MSE, and covariance structures, to an earlier multivariate estimator defined in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF]Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] and also to the univariate estimator defined in [START_REF] Wendt | [END_REF]. The potential benefits of the proposed eigen-wavelet estimation procedures are illustrated by means of epileptic seizure prediction from multichannel EEG data. A toolbox and its documentation for both the synthesis of multivariate selfsimilar processes and for the estimation of the vector of selfsimilarity parameters will be made publicly available at the time of publication.

II. MULTIVARIATE SELFSIMILARITY: STATE OF THE ART

A. Model: Operator Fractional Brownian Motion (ofBm) 1) Definition: Operator Fractional Brownian Motion, originally proposed in [START_REF] Mason | Sample path properties of similar Gaussian fields[END_REF]- [START_REF]Exponents, symmetry groups and classification of operator frac-Brownian motions[END_REF], is defined as a jointly Gaussian M-variate process, via an harmonizable representation,

B H,A (t) := R e itf -1 if f -H+ 1 2 I + A + f -H+ 1 2 I - A B(df ), (1)
where I denotes the M × M identity matrix, B(df ) a complex-valued Gaussian Hermitian random measure, such that E B(df ) B(df ) = df (with the Hermitian transposition), and A and H two M × M complex-valued matrices.

Furthermore, hereinafter it is also assumed that [START_REF] Pipiras | Integral representations and of fractional Brownian motions[END_REF]: ASSUMPTION OFBM1: The Hurst matrix has the form

H = Wdiag(H 1 , . . . , H M )W -1 , (2)
with W a M × M real-valued invertible matrix, and where

0 < H 1 ≤ . . . ≤ H M < 1. (3)
ASSUMPTION OFBM2: AA is positive definite, i.e., det (AA ) > 0. ASSUMPTION OFBM3: AA is a real-valued matrix. OFBM2 ensures that B H,A is a well-defined process ; OFBM3 restricts B H,A to the special case of time reversible (or covariant under the change t → -t) statistical properties ; W in OFBM1 is referred to as the mixing matrix and provides an important source of versatility of the model.

From Definition (1), the covariance matrix is obtained as:

E B H,A (t)B H,A (s) = R (e itf -1)(e -isf -1) |if | 2 f -(H-1 2 I) + AA f -(H -1 2 I) + +f -(H-1 2 I) - AA f -(H -1 2 I) - df, (4)
where E stands for the ensemble average.

2) Multivariate selfsimilarity: Importantly, ofBm satisfies the general multivariate selfsimilarity relation [START_REF] Pipiras | Integral representations and of fractional Brownian motions[END_REF] [START_REF] Palmer | Fractal geometry: a tool for describing spatial patterns of plant communities[END_REF] where a H = k≥0 log k a H k /k! and f dd = denotes equality for all finite-dimensional distributions. In essence, Eq. ( 5) states that the statistics of B H,A (t) are covariant under a change of time scale, t → t/a, up to a change in amplitude B → a H B. Importantly, this holds for any dilation factor a > 0. However, the factor controlling the change in amplitude → a H involves jointly all entries of the vector of selfsimilarity parameters H = (H 1 , . . . , H M ), through the mixing matrix W.

{B H,A (t)} t∈R f dd = {a H B H,A (t/a)} t∈R , ∀a > 0,
Of particular pedagogical interest is the case when the scaling exponent matrix H is diagonal, H = diag(H 1 , . . . , H M ) (equivalently, W = I): M -variate selfsimilarity in Eq. ( 5) then simplifies into a collection of M , component-wise, univariate selfsimilarity relations, each controlled by a single selfsimilarity parameter, i.e., ∀m = 1, . . . , M and a > 0:

{(B H,A ) m (t)} t∈R f dd = {a Hm (B H,A ) m (t/a)} t∈R . (6)
B. Eigen-wavelet analysis 1) Multivariate wavelet spectrum: Let ψ 0 denote the socalled mother wavelet defined as a compact support, unit norm reference pattern R ψ 0 (t) 2 dt = 1, characterized by its number of vanishing moments N ψ , and satisfying regularity conditions quantified by the decay of its Fourier transform ψ0 : There is α > 1 such that sup

f ∈R ψ0 (f ) (1 + |f |) α < +∞.
Let Y denote a M × N multivariate process, with M and N the numbers of components and samples. The multivariate Discrete Wavelet Transform coefficients at scale

2 j consist of a M × N matrix D Y (j, k) = R ψ(j, k)Y (t)dt, with ψ(j, k) = 2 -j/2 ψ 0 (2 -j t -k) a collection of
dilated and translated templates of ψ 0 . For further details on wavelet transforms, the reader is referred, e.g., to [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

The wavelet spectrum is defined as the collection of M ×M scale-dependent random matrices:

S Y (2 j ) = 1 n j nj k=1 D Y (2 j , k)D Y (2 j , k) . (7)
2) Multivariate wavelet spectrum for ofBm: For B H,A , it has been shown that wavelet coefficients reproduce selfsimilarity as [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF]Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]:

D B H,A (2 j , k) k∈Z f dd = 2 j(H+ 1 2 I) D B H,A (2 0 , k) k∈Z . (8)
The wavelet spectrum thus becomes:

S B H,A (2 j ) d = 2 j(H+ 1 2 I) S B H,A (2 0 ) 2 j(H + 1 2 I) . (9)
When W = I, Eq. [START_REF] He | Scale-free brain activity: past, present, and future[END_REF] shows that (the expectation of) each entry of S B H,A (2 j ) generally consists of a mixture of power laws with respect to the scales 2 j . In other words,

ES m,m (2 j ) = M k=1 M k =1 α m,m k,k 2 j(H k +H k +1) , (10)
with α m,m k,k a priori depending on both H and A. 3) Estimation of the vector of selfsimilarity parameters: It has been shown in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF]Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] that each eigenvalue λ m (2 j ) of ES B H,A (2 j ) behaves asymptotically as a power law with respect to scales 2 j , with scaling exponent controlled by H m ,

λ m (2 j ) ∼ ξ m (2 0 ) (2 j ) 2Hm+1 , as j → +∞. (11)
Eq. ( 11) naturally suggests that the estimation of the vector of selfsimilarity parameters H can be performed by a collection of individual, for each m = 1, . . . , M , linear regressions on the logarithms of the eigenvalues λm (2 j ) of S B H,A (2 j ) versus the logartihms of the scales 2 j :

ĤM m = 1 2   j2 j=j1 w j log 2 λm (2 j ) -1   , (12)
where j ∈ (j 1 , j 2 ) denote the range of octaves involved and the regression weights satisfy: j2 j=j1 jw j = 1 and j2 j=j1 w j = 0. The asymptotic estimation performance was thoroughly studied theoretically in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF]Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF].

III. MULTIVARIATE FRACTIONAL BROWNIAN MOTION A. Definition

For the specific nonmixing case (W = I), the selfsimilarity relation provided in Eq. ( 6) suggests, as a first contribution of this work, the construction of another type of multivariate selfsimilar process. It is hereinafter referred to as multivariate fractional Brownian motion (mfBm) B Σ,W,H (t) and defined as a mixture of correlated univariate fractional Brownian motions:

B Σ,W,H (t) := WB Σ,H (t) , (13)
where W denotes a M ×M real-valued invertible matrix. Also, From the definition of univariate fBm [START_REF] Samorodnitsky | Stable non-Gaussian Random Processes: Stochastic Models with Infinite Variance[END_REF], we can show that a mfBm B Σ,W,H is a zero-mean multivariate Gaussian process with covariance function given by: 4) and ( 14), calculations not reported here yield the relations:

B Σ,H (t) = [B Σ,H1 (
E B Σ,W,H (t) B Σ,W,H (s) = 1 2 W |t| diag(H) Σ |t| diag(H) +|s| diag(H) Σ |s| diag(H) -|t -s| diag(H) Σ |t -s| diag(H) W . ( 14
)
AA * = W • (G Σ) • W and H = W diag(H) W -1 . (15)
In [START_REF] Abry | Multiscale Anisotropic Texture Analysis and Classification of Photographic Prints: Art scholarship meets image processing algorithms[END_REF], G is a M × M matrix with entries g m,m = Γ(H m + H m + 1) sin((H m + H m )π/2)/(2π), where Γ is the standard Gamma-Euler function and the symbol stands for the Hadamard matrix product. Eq. [START_REF] Abry | Multiscale Anisotropic Texture Analysis and Classification of Photographic Prints: Art scholarship meets image processing algorithms[END_REF] shows that the mfBm B Σ,W,H is a special case of the ofBm B H,A , satisfying Assumptions OFBM1 to OFBM3. The mfBm B Σ,W,H thus provides a model for multivariate selfsimilarity. It can be considered a practical counterpart of the formal ofBm model, yet better suited for signal processing purposes and applications to real-world data. It splits the abstract matrices A and H from the ofBm model into interpretable ingredients that are closer to applications: a vector of scaling exponents H, a (premixing or intrinsic) covariance matrix Σ, and a mixing matrix W.

In particular, Eq. [START_REF] Abry | Multiscale Anisotropic Texture Analysis and Classification of Photographic Prints: Art scholarship meets image processing algorithms[END_REF] shows that mfBm B Σ,W,H necessarily fulfills Assumption OFBM3 and thus implies by construction time-reversible statistics. An alternative permitting non time reversible statistics had been considered in [20], yet not allowing for mixing hence less versatile for applications. Furthermore, the mfBm model allows us to better grasp a key feature of multivariate selfsimilarity. Starting with the simple bivariate setting, M = 2, combining Assumption OFBM2 and Eq. ( 15) yields:

ρ 2 12 ≤ ρ 2 max := Γ(2H 1 + 1)Γ(2H 2 + 1) sin(πH 1 ) sin(πH 2 ) (Γ(H 1 + H 2 + 1) sin(π/2(H 1 + H 2 ))) 2 . ( 16
) This shows that the correlation between components and differences in selfsimilarity parameters cannot be chosen independently. This is illustrated in Fig. 1: The largest allowed (squared) cross-correlation amongst components decreases when the difference |H 2 -H 1 | increases, and conversely. These calculations extend to M-variate settings with larger (absolute values of) correlations amongst components restricting the largest achievable discrepancies |H m -H m | amongst entries of the vector of selfsimilarity parameters (see also [20]).

Finally, B Σ,H ≡ B Σ,W=I,H is obviously a special case of mfBm (cf., e.g., [START_REF] Wendt | [END_REF]- [START_REF] Coeurjolly | Wavelet the multivariate Brownian motion[END_REF]). Yet, it displays no mixing across components, rendering the model less versatile for practical purposes.

C. Wavelet analysis of mfBm

For mfBm, Eqs. ( 8) and ( 9) take the specific forms

{D B Σ,W,H (2 j , k)} k∈Z f dd = W 2 j(H+ 1 2 ) D B Σ,H (2 0 , k) k∈Z , (17) S B Σ,W,H (2 j ) d = W 2 j(H+ 1 2 ) S B Σ,H (2 0 ) 2 j(H+ 1 2 ) W , (18)
with 2 j(H+ 1 2 ) = diag(2 j(H1+1/2) , . . . , 2 j(H M +1/2) ). This disentangles the contributions of the dilation vs. mixing operators to the shaping of the wavelet coefficients and spectrum.

Furthermore, Eqs. ( 10) and ( 11) obviously hold for mfBm. This permits to show, respectively, that both α m,m k,k and ξ m (2 0 ) depend on W and Σ, but not on H. An extension of the calculations in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF]Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] shows that ξ m (2 0 ) actually depends only on the entries

S B Σ,W,H (2 0 ) (m ,m ) for m ≤ m , m ≤ M (see Remark A.1 in Section A of the supplementary material).
Finally and importantly, Eq. ( 11) demonstrates that the wavelet eigenvalue representations for both mfBm and ofBm disentangle the mixed power laws appearing in the wavelet spectrum S B Σ,W,H (2 j ). Thus, they restore the close relationship between power laws and scale-free dynamics, which had gotten altered by the mixing effect of the matrix W. This motivates the estimation procedure ĤM in Eq. ( 12), which thus applies to mfBm as well as to ofBm.

IV. SELFSIMILARITY PARAMETER VECTOR ESTIMATION

A. Repulsion bias-corrected multivariate estimation 1) Repulsion bias:

ĤM in Eq. ( 12) is based on eigenvalues computed from the wavelet spectra S B H,A (2 j ), namely, on the empirical wavelet coefficient covariance matrices estimated at different scales. However, it is well known that, for finitesize sample estimation of covariance matrices, eigenvalues undergo the so-called repulsion effect: Estimated eigenvalues tend to be farther apart one from the others, compared to true eigenvalues, all the more as the sample size is decreasing compared to the number of components (cf. e.g., [START_REF] Tao | Topics in Random Matrix Theory[END_REF]).

In the wavelet spectrum representation, this repulsion effect is enhanced by the fact that empirical wavelet coefficient covariance matrices are estimated with sample size n j that depends on analysis scales 2 j , essentially as n j = 2 -j N . This induces scale-dependent repulsion biases in λm (2 j ), and hence an overall bias in ĤM [START_REF] Lucas | Bootstrap for testing the equality of selfsimilarity exponents across multivariate time series[END_REF]. The repulsion bias is illustrated in Fig. 2 (left), which displays log 2 λm (2 j ) vs log 2 2 j = j. It is based on a synthetic mfBm with M = 6 components of size N = 2 12 , with H 1 = . . . = H 6 = 0.6, a randomly selected orthonormal mixing matrix W = I and a correlation matrix Σ with off-diagonal coefficients set to 0.8. The five smallest curves should thus theoretically superimpose, whereas discrepancies of amplitudes in estimated eigenvalues, increasing with scales, can be observed.

2) Repulsion bias correction: Elaborating on [START_REF] Lucas | Bootstrap for testing the equality of selfsimilarity exponents across multivariate time series[END_REF], a novel estimation procedure for the estimation of H is devised and studied, hereinafter referred to as the bias-corrected estimator ĤM,bc = ( ĤM,bc 1 , . . . , ĤM,bc

M

). The main idea behind its construction lies in using the same number of wavelet coefficients in the wavelet spectrum estimation at each scale. As a result, the eigenvalues computed from the wavelet spectra at different scales 2 j undergo an equivalent strength repulsion effect, i.e., one that does not vary with the scale 2 j . To that end, let 2 j2 denote the largest analysis scale used in estimation (cf. Eq. ( 12)) and let n j2 be the number of corresponding wavelet coefficients at that scale. At each scale 2 j1 ≤ 2 j ≤ 2 j2 , a collection of wavelet spectra are computed from τ = 1, . . . , 2 j2-j consecutive non-overlapping windows of wavelet coefficients with same size n j2 :

S (τ ) (2 j ) := 1 n j2 τ nj 2 k=1+(τ -1)nj 2 D B H,A (2 j , k)D B H,A (2 j , k) T .
(19) At each scale, the logarithms of the eigenvalues λ(τ) m (2 j ) computed from each wavelet spectrum S (τ ) (2 j ) are then averaged across windows:

log 2 λ m (2 j ) := 2 j-j2 2 j 2 -j τ =1 log 2 λ(τ) m (2 j ). ( 20
)
The impact of this bias correction procedure on the scaling of the eigenvalues is illustrated in Fig. 2 (right), which clearly shows that the theoretical slopes 2H m + 1 are much better reproduced by the log 2 λ m (2 j ) than by the classical log 2 λm (2 j ). The bias-corrected multivariate estimates ĤM,bc m are thus defined as linear regressions across (the logarithms of the) scales 2 j of the averaged log-eigenvalues log 2 λ m (2 j ):

ĤM,bc m = 1 2   j2 j=j1 w j log 2 λ m (2 j ) -1   . (21)
The asymptotic estimation performance is studied theoretically in Section IV-B. The finite-size estimation performance is investigated numerically in Section V.

3) Nonmixing case: In the nonmixing case, i.e., when W = I, the selfsimilarity matrix reduces to H = diag(H).

While the multivariate estimators ĤM m in Eq. ( 12) and ĤM,bc m in Eq. ( 21) obviously apply to the nonmixing case, a natural estimation procedure, extending univariate analysis, consists in performing linear regressions on the logarithms of the diagonal coefficients S m,m (2 j ) only:

ĤU m = 1 2   j2 j=j1 w j log 2 S m,m (2 j ) -1   . ( 22
)
The resulting estimators ĤU m , intrinsically univariate in spirit, were intensively studied in [START_REF] Wendt | [END_REF] for nonmixing mfBm.

B. Theoretical study of asymptotic estimation performance 1) Asymptotic framework: The estimation performance is studied theoretically in the asymptotic double and joint limits of large sample sizes N → +∞, with linear regressions performed at coarse scales, in the range (j 1 (N ), j 2 (N )) → +∞. Technically, the scaling range is defined as

(j 1 (N ) = j 0 1 +log 2 a(N ), j 2 (N ) = j 0 2 +log 2 a(N ))
, where (j 0 1 , j 0 2 ) is an arbitrarily chosen range for N 0 , and where a(N ) must satisfy:

a(N )2 j2 ≤ N, a(N ) N + N a(N ) 2 +1 -→ 0, as N → ∞, (23)
= min min

{2≤i≤M | Hi-Hi-1>0} (H i -H i-1 ), H 1 2 + 1 4 .
(24) This essentially implies that log 2 a(N ) increases as β log 2 N , with 1/(2 + 1) < β < 1. This also implies that the number of scales involved in linear regressions does not vary with sample size N : j 2 (N ) -j 1 (N ) + 1 = j 0 2 -j 0 1 + 1. 2) Asymptotic consistency: Theorem IV.1 establishes the asymptotic consistency of the bias-corrected estimator ĤM,bc :

Theorem IV.1. Under Assumptions (OFBM1-3) and for a mother wavelet as defined in Section II-B1, (with P → the convergence in probability):

(i) For any j = j 0 1 , . . . , j 0 2 , τ = 1, . . . , 2 j 0 2 -j and m = 1, . . . , M , as N → ∞,

λ(τ) m (a(N )2 j ) a(N ) 2Hm+1 P → ξ m (2 j ) = (2 j ) 2Hm+1 ξ m (2 0 ) > 0, (25) where λ(τ) 1 (a(N )2 j ), . . . , λ(τ) M (a(N )2 j ) are the eigen- values of S (τ ) (a(N )2 j ). Likewise, as N → ∞, λ (τ ) m (a(N )2 j ) a(N ) 2Hm+1 → ξ m (2 j ) = (2 j ) 2Hm+1 ξ m (2 0 ) > 0, (26) where λ
(τ ) 1 (a(N )2 j ), . . . , λ (τ )
M (a(N )2 j ) are the eigenvalues of ES (τ ) (a(N )2 j ).

(ii) Let a(N ) be as in [START_REF] Pipiras | Integral representations and of fractional Brownian motions[END_REF]. Then, as N → ∞,

ĤM,bc

m P → H m , ∀m = 1, . . . , M, (27)
Proof. For (i), the proof follows from adapting the technique for establishing Theorem 3.1 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], see Section A of the supplementary material for details.

For (ii), Eq. ( 25) leads to the following equivalence:

ĤM,bc m = 1 2 j 0 2 j=j 0 1 w j 2 j 0 2 -j 2 j 0 2 -j τ =1 log 2 λ(τ) m (a(N )2 j ) a(N ) 2Hm+1 (28)
-log 2 a(N ) 2Hm+1 -1 2

P ∼ 1 2 j 0 2 j=j 0 1 w j 2 j 0 2 -j 2 j 0 2 -j τ =1 log 2 (2 j ) 2Hm+1 ξ m (2 0 ) -log 2 a(N ) 2Hm+1 - 1 2 P ∼ 1 2 j 0 2 j=j 0 1 w j (2H m + 1)j + log 2 ξ m (2 0 ) a(N ) 2Hm+1 - 1 2
.

Hence, by the definition of w j , as N → ∞, ĤM,bc

m P ∼ 1 2 [(2H m + 1)] - 1 2 = H m . (29)
This establishes [START_REF] Lucas | Bootstrap for testing the equality of selfsimilarity exponents across multivariate time series[END_REF].

3) Simple eigenvalues condition: Our theoretical study of the asymptotic distribution relies on the condition that, asymptotically, all wavelet spectra eigenvalues are simple. CONDITION (C0):

∀m 1 , m 2 ∈ {1, . . . , M }, with m 1 = m 2 , H m1 = H m2 ⇒ ξ m1 (2 0 ) = ξ m2 (2 0 ), (30)
which states that either the selfsimilarity parameters H m are different, or when they are equal, H m1 = H m2 , m 1 = m 2 , the constants ξ m (2 0 ) in Eq. ( 11) are different. 4) Asymptotic normality: Theorem IV.2 establishes the asymptotic normality of the bias-corrected estimator ĤM,bc :

Theorem IV.2. Let n a,j = N a(N )2 j , j = j 0 1 , . . . , j 0 2 , be the number of wavelet coefficients available at scale a(N )2 j . Under Assumptions (OFBM1-3), for a mother wavelet as defined in Section II-B1, and assuming that (C0) holds, with a(N ) defined in [START_REF] Pipiras | Integral representations and of fractional Brownian motions[END_REF] and where d → denotes the convergence in distribution:

(i) Then, as N → ∞, n a,j 0 2 log 2 λ(τ) m (a(N )2 j ) (31)
-

log 2 λ (τ ) m (a(N )2 j ) j=j 0 1 ,...,j 0 2 m=1,...,M τ =1,...,2 (j 0 2 -j 0 1 )-j d → N (0, Υ B ),
where Υ B is a symmetric positive semidefinite matrix.

(ii) Moreover, as N → ∞, N/a(N ) ĤM,bc m -H m m=1,...,M d → N (0, Σ B ), (32)
for some Σ B ∈ S ≥0 (M, R).

Proof. The proof follows from adapting the proof of asymptotic normality for the non bias-corrected multivariate estimator ĤM from Theorem 3.2 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], as detailed in Section B of the supplementary material.

V. EMPIRICAL FINITE-SIZE ESTIMATION PERFORMANCE

The goal of the present section is to complement the theoretical study of estimation performance in the limit of large scales and large sample sizes. We do so by investigating their finite sample size counterparts by means of Monte Carlo simulations, conducted on synthetic mfBm. Different parametric configurations (Σ, W, H) are used so as to further assess their impact on estimation performance.

The performance of the proposed bias-corrected multivariate estimator ĤM,bc (Eq. ( 21)) is compared, in terms of biases, variances, Mean Squared Error (MSE), and covariance structures, to those of the multivariate estimator ĤM (Eq. ( 12)) and of the univariate estimator ĤU (Eq. ( 22)).

A. Numerical simulation set-up

The Monte Carlo simulations involve N MC = 1000 independent realizations of synthetic mfBm. The increments {B Σ,W,H (t+1)-B Σ,W,H (t)} t=1,...,N of a mfBm, referred to as multivariate fractional Gaussian noise, form a multivariate stationary Gaussian process with covariance directly stemming from Eq. ( 14). Thus, it can be synthesized numerically by circulant embedding, cf. e.g., [START_REF] Helgason | Fast and exact synthesis of stationary multivariate gaussian time series using circulant embedding[END_REF].

Results reported here are restricted to M = 6 for the sake of clarity ; additional results for M = 20, available upon request, led to identical conclusions. The range of sample sizes covered is N ∈ {2 13 , . . . , 2 18 }. The parameters Σ, W and H are chosen to investigate a set of representative configurations: A generic covariance matrix Σ is chosen: Σ m,m = 0.7 |m-m | , 1 ≤ m ≤ m ≤ M . A mixing matrix W is randomly drawn from the set of M × M real-valued nonorthogonal matrices and kept fixed for all experiments. The selfsimilarity parameter vector H = (H 1 , . . . , H 6 ) is either set to H m = 0.7, ∀m, or H = (0.4, 0.5, 0.6, 0.65, 0.7, 0.8). These different parameter configurations are chosen so that CONDI-TION (C0) holds -this condition is required for asymptotic normality claims (cf. Section IV-B3). Fig. 3: Finite-size joint normality for ĤM,bc . Quantilequantile plots of the empirical distribution of the Mahalanobis distance computed from the bias-corrected multivariate estimates ĤM,bc against a χ 2 distribution with M degrees of freedom (expected theoretically under joint normality), for different sample sizes N and different configurations.

Config1 (Σ = I, W = I, H m all different)
Wavelet analysis is performed using the least asymmetric Daubechies2 mother wavelet, with N ψ = 2 vanishing moments [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Linear regressions are performed, with w j as defined in [START_REF] Wendt | [END_REF], over scales that depend on the sample size N , as analyzed theoretically in Section IV-B1. Namely, we pick 2 j1 = a(N )2 6 to 2 j2 = a(N )2 9 , where a(N ) = 2 β log 2 (N/N0) with β = 0.9 and N 0 = 2 13 (cf. Eq. ( 23)).

B. Asymptotic normality assessment

The normality of ĤM,bc is assessed by computing quantilequantile plots for the squared Mahalanobis distance statistic

( ĤM,bc -E ĤM,bc )Var( ĤM,bc ) -1 ( ĤM,bc -E ĤM,bc ) . (33)
In [START_REF] Lucas | Counting the number of different scaling exponents in multivariate scale-free dynamics: Clustering by bootstrap in the wavelet domain[END_REF], the ensemble average E is practically computed as the mean across independent realizations. Samples of the random variable [START_REF] Lucas | Counting the number of different scaling exponents in multivariate scale-free dynamics: Clustering by bootstrap in the wavelet domain[END_REF] are plotted against a χ 2 distribution with M degrees of freedom, since this distribution is expected to hold under the exact joint normality of ĤM,bc . Fig. 3 reports such quantile-quantile plots, for the different configurations and for different sample sizes. It shows very satisfactory straight lines, even at the smallest sample size N = 2 13 . This suggests that the joint normality for ĤM,bc holds very satisfactorily for finite sample size estimation, even for small sample sizes. It thus indicates fast convergence to asymptotic normality, proven theoretically in Theorem IV.2 (Eq. ( 32)).

C. Estimation performance: bias, covariance, MSE

For any estimate of the vector of selfsimilarity parameters Ĥ, matrices of bias, covariance and MSE are defined as The estimation performance is quantified by the spectral norm of these matrices, i.e., by the largest absolute value of their eigenvalues [START_REF] Meyer | Matrix Analysis and Applied Linear Algebra[END_REF].

Bias 2 ( Ĥ) = (E Ĥ -H)(E Ĥ -H) , (34)
Cov( Ĥ) = E (E Ĥ -Ĥ)(E Ĥ -Ĥ) , (35)
MSE( Ĥ) = E ( Ĥ -H)( Ĥ -H) . (36)
Fig. 4 reports the (logarithms of the) spectral norms of the bias, covariance and MSE for the proposed bias-corrected multivariate estimator ĤM,bc , compared to those for the classical multivariate estimator ĤM and for the univariate estimator ĤU . The results are shown as functions of the logarithm of the sample size N and for each of the different configurations.

For the generic case represented by Config1 (Fig. 4, first row), the univariate estimator ĤU suffers from significant bias stemming from the component mixing caused by W. Thus, it is significantly outperformed by the multivariate estimators ĤM and ĤM,bc . The latter show equivalent performance as the combination of Σ, W and different H m yield different eigenvalues for ES(2 j ) at all scales and thus small repulsion bias in their estimated values.

For Config2 (Fig. 4, second row), with mixing yet equal H m , the eigenvalues of ES(2 j ) lie much closer to each other at most scales and hence suffer from substantial repulsion biases in their estimation. This results in a significant bias for the classical multivariate estimator ĤM . Fig. 4, middle row, shows that the proposed bias-corrected multivariate estimator ĤM,bc succeeds in significantly reducing bias without an increase in the variance. This leads to an overall much improved estimation performance as quantified by the MSE. While the univariate estimator ĤU displays small bias, as expected in view of the equality of all selfsimilarity parameters, it also shows greater variance, thus, greater MSE, when compared to the proposed ĤM,bc . For the nonmixing case, Config3 (Fig. 4, third row), the univariate estimator ĤU naturally shows much smaller biases compared to the multivariate estimators ĤM and ĤM,bc . Yet, it also displays larger variances. Thus, as a result of bias-variance trade-off, the MSE of all three estimators are equivalent. Importantly, this shows that even for data corresponding to nonmixing situations, there is no cost in estimation performance associated with the use of multivariate estimators.

For the nonmixing case with equal H m , Config4 (Fig. 4, fourth row), the univariate estimator ĤU and bias-corrected multivariate estimator ĤM,bc perform similarly to Config2. This suggests that, for H m all equal, the presence of mixing has no impact in terms of bias or variance. However, estimation performance of the multivariate estimator ĤM is slightly more affected than in Config2, showing that the repulsion bias is more substantial in the absence of mixing.

Finally, for the mixing case with different H m and no correlation, Config5 (Fig. 4, fifth row), the performances are similar to those for Config1 for all three estimators, indicating their robustness against component correlations Σ.

All together, the results reported in Fig. 4 lead us to conclude that the proposed bias-corrected multivariate estimator ĤM,bc displays the best estimation performance regardless of the actual configuration of the parameters. In practice it is not a priori known whether or not mixing is present, or whether selfsimilarity parameters are pairwise distinct or equal. We have seen that the estimator ĤU fails when mixing is present. Also, we have found out that the estimator ĤM fails when repulsion bias is present, i.e., when one needs to estimate identical selfsimilarity parameters. By contrast, ĤM,bc displays neither one of these limitations -this versatility is clearly an important property that motivates its practical use.

D. Covariance structure of

ĤM,bc 1) Correlation structure: We complement the estimation performance assessment by turning to the dependence structure amongst estimates -the latter constitutes a key feature in the practical use of the estimators of the selfsimilarity parameters. Asymptotic normality, as well as finite sample size approximate normality, hint on the importance of the correlation structure amongst estimated selfsimilarity parameters. compared to the multivariate estimators (cf. Fig. 4 above). These empirical observations are complemented and framed by the analytical calculations detailed in Sections C and D of the supplementary material. These show the following striking fact. Suppose CONDITION (C0) holds, and further assume that the weak correlation structure amongst wavelet coefficients of mfBm can be approximated into exact decorrelation. Then, the covariance matrices Σ B in Eq. ( 31) and Υ B in Eq. ( 32) of Theorem IV.2 are well approximated, to the first order, by diagonal matrices. These empirical observations and analytical calculations reveal an important advantage of the mutivariate ĤM,bc (and ĤM ) against the univariate ĤU , even for nonmixing data: namely, that the former consist of asymptotically Gaussian, weakly dependent random vectors. This is a major practical feature, e.g., in the design of tests for the equality of selfsimilarity parameters, which is of significant interest in many applications [START_REF] Lucas | Bootstrap for testing the equality of selfsimilarity exponents across multivariate time series[END_REF], [START_REF] Lucas | Counting the number of different scaling exponents in multivariate scale-free dynamics: Clustering by bootstrap in the wavelet domain[END_REF].

2) Variance: Besides approximate asymptotic decorrelation, the diagonal entries of the matrices Σ B and Υ B remain to be discussed. The calculations reported in Remark A.2 in Section B, and in Sections C and D of the supplementary material show that these diagonal entries can be first-order approximated by a constant:

Υ B 2 (log 2 e) 2 I and Σ B (log 2 e) 2 2 j 0 2 j=j 0 1 w 2 j 2 j I. (37)
Specifically, this implies that, under the approximation of decorrelation between wavelet coefficients of mfBm and under (C0), the finite-size variances of ĤM,bc m and ĤM m , for every m ∈ {1, . . . , M }, can be approximated as: ) decrease with sample size N and depend on the range of regressions scales j 1 and j 2 . However, to the first order, these depend on none of the model parameters Σ, W and H. This renders them comparable for all components. Interestingly, it was shown in [START_REF] Wendt | [END_REF] that, for the univariate ĤU m in nonmixing settings, the variance can be approximated to the first order by the same closed-form relation: Var( ĤU m ) ≈ V N . To further assess the practical relevance of the finite sample size approximation in Eq. 38, Fig. 6 reports the relative difference between the empirical variances of ĤU m , ĤM m and ĤM,bc m and V N as a function of the logarithm of the sample size N for m = 1, . . . , 6 and for the different configurations. Fig. 6 (left column), first, shows that, in all configurations, the variance of the univariate estimator ĤU m is approximated by V N with a 20% relative error resulting from the non-exact decorrelation of eigenvalues between scales. This first result extends the results reported in [START_REF] Wendt | [END_REF], limited to nonmixing cases. Fig. 6, second, confirms that the variances of the three estimators neither depend on H (Config1 vs. Config2), nor on W (Config1 vs. Config3), nor on Σ (Config1 vs. Config5). Fig. 6, third, confirms that the variances of the multivariate estimators are, on average, lower than those of the univariate estimator. In summary, Fig. 6 shows that the closed-form approximation V N in Eq. ( 38) is valid in all configurations for the three estimators, even for small sample sizes. This constitutes a valuable result for practical use in applications e.g., when testing the equality of selfsimilarity parameters.

Var( ĤM,bc m ) ≈ Var( ĤM m ) ≈ (log 2 e) 2 2 j2 j=j1 w 2 j n j := V N . (38)

VI. APPLICATION TO EPILEPTIC SEIZURE PREDICTION

Selfsimilarity has been broadly used to model temporal dynamics in macroscopic brain activity. The purpose of the modeling effort is either to study healthy brain mechanisms (cf. e.g., [START_REF] He | Scale-free brain activity: past, present, and future[END_REF], [START_REF] La Rocca | Self-similarity and multifractality in human brain activity: a waveletbased analysis of scale-free brain dynamics[END_REF], [START_REF] Ciuciu | Interplay between functional connectivity and scale-free dynamics in intrinsic fMRI networks[END_REF]), or to diagnose brain disorders and pathologies. Notably, there is a particular interest in epileptic seizure prediction [START_REF] Domingues | Multifractal analysis for cumulant-based epileptic seizure detection in eeg time series[END_REF]. Yet, most studies have so far remained based on univariate analysis of scale-free dynamics, while brain dynamics is naturally conducive to multivariate analysis. The aim of this section is, thus, to comparatively assess the potential of multivariate and univariate selfsimilarity parameter estimations in the detection of preictal states (i.e., time periods prior to seizure onsets) against interictal states (i.e., periods distant in time from any seizure).

The multi-channel EEG recordings studied here were collected at the Boston Children's Hospital from pediatric patients with medically intractable seizures, from the CHB-MIT Scalp EEG database [START_REF] Goldberger | PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals[END_REF], available at https://physionet.org/content/ chbmit/1.0.0/. EEG recordings stem from the International 10-20 system of EEG electrode positions and nomenclature, last for at least one hour, and are sampled at 256Hz. Multivariate analyses are performed for the 19 channels, within sliding 2minute long windows. Wavelet analysis is performed using the least asymmetric Daubechies2 mother wavelet. Fig. 7 reports, for both interictal and preictal windows, the logarithms of the diagonal entries of the wavelet spectrum S(2 j ), log 2 S m,m (2 j ), the log-eigenvalues log 2 λm (2 j ) and the proposed bias-corrected log-eigenvalues log 2 λ m (2 j ), all of them as functions of scales 2 j for m = 1, . . . , 19. Fig. 7 shows linear behaviors for these functions across fine scales 1 ≤ 2 j ≤ 2 7 , that hence indicate scale-free dynamics in these EEG time series, at frequencies above 1Hz, in agreement with the literature on epilepsy, cf. e.g., [START_REF] Gadhoumi | Scale invariance properties of intracerebral EEG improve seizure prediction in mesial temporal lobe epilepsy[END_REF]. These scalefree dynamics hold for both preictal and interictal states, and both component-wise (univariate, log 2 S m,m (2 j )) and jointly (multivariate, log 2 λm (2 j ) and log 2 λ m (2 j )). The observed multivariate selfsimilarity constitutes per se a significant outcome of this work as selfsimilarity in EEG time series had so far only been observed component per component and not jointly. Combined with literature suggesting that differences between preictal and interictal dynamics are related to frequencies ranging from 10Hz to 85Hz [START_REF] Gadhoumi | Scale invariance properties of intracerebral EEG improve seizure prediction in mesial temporal lobe epilepsy[END_REF], these observations suggest performing linear regressions for the estimation of H across scales ranging from 2 1 ≤ 2 j ≤ 2 4 . Fig. 8 compares the (logarithms of) sorted Wilcoxon ranksum test p-values p τ (m) measuring the significance of the differences between the medians of preictal vs. interictal distributions, for the univariate estimator ĤU (Eq. ( 22)), the multivariate estimator ĤM (Eq. ( 12)) and the bias-corrected multivariate estimator ĤM,bc (Eq. ( 21)). The p-values are ranked and compared to Benjamini-Hochberg multiple hypothesis thresholds [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF], computed at a false discovery rate α = 0.05. For all three estimators, the medians of the distributions of estimated selfsimilarity parameters clearly differ from preictal to interictal states, thus indicating that the scale-free structure of the temporal dynamics is modified by forthcoming epileptic seizures. Moreover, the bias-corrected multivariate Fig. 7 importance, the distributions of selfsimilarity parameters are found to differ from to another. This demands a procedure to detect preictal vs. interictal states. Interested readers are referred to the in-depth study carried out [START_REF] Lucas | Epileptic seizure prediction from eigen-wavelet multivariate selfsimilarity analysis of multi-channel EEG signals[END_REF].

VII. CONCLUSION

In this work, we have, first, proposed a multivariate selfsimilarity model, defined as a linear mixture of correlated fractional Brownian motions, well-suited to adjust the versatility of real-world multivariate data, and formally related it to a former mathematical model. Second, we have turned the formal proposition of wavelet spectrum eigenvalue-based estimation for the vector of selfsimilarity parameters proposed in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF]Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF], into a versatile and robust bias-corrected estimation procedure that can actually be used on real-world data. Its performance was studied mathematically in terms of asymptotics and by numerical simulations for finite-size samples. This permitted showing that the proposed multivariate estimators work in a broader context (mixing) than the univariate one and still performs better (decorrelation) in the nonmixing case, where the univariate estimator can be used. Third, its potential was illustrated at work on real multi-channel EEG data in a context of epileptic seizure prediction. A documented toolbox implementing estimation will be made publicly available at the time of publication.

Future investigations include incorporating a multivariate time-scale bootstrap procedure so as to assess confidence levels in estimation from a single realization of data. Also, largedimensional asymptotics will be studied both theoretically and numerically to account for situations where the number of components grows as fast as sample size, of major practical importance in real-world applications.

where u m is a unit eigenvector of M associated with λ m (M) (see [START_REF] Magnus | On differentiating eigenvalues and eigenvectors[END_REF], p. 182, Theorem 1). Now, for every B ∈ S ≥0 (M, R), define the function f N,m by f N,m (B) = log 2 λ m Wa diag(H) Ba diag(H) W * a 2Hm . (51)

For any fixed δ > 0, let O δ = {B ∈ S ≥0 (M, R) : B-B(2 j ) < δ}. By Proposition 3.1 in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], det B(2 j ) > 0. Thus, an application of the proof of Theorem 3.2 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF] implies that there exists some δ 0 > 0 such that, for large enough N and for any B ∈ O δ0 , the mean value theoremtype relation x * ,m,i1 •x * ,m,i2 • π i1,i2 B (τ ) a (2 j ) -B(2 j ) .

(58) Now note that (58) further holds for all τ , m and j, and also that these expansions start from the same sequence of asymptotically Gaussian random matrices B (τ ) a (2 j ) (see (48)). Hence, the jointly asymptotically normal limit (31) holds. Thus, (i) is proven.

Statement [START_REF] Meyer | Matrix Analysis and Applied Linear Algebra[END_REF] can be proven by naturally adapting the proof of Theorem 3.2 in [START_REF] Abry | Wavelet eigenvalue regression in high dimensions[END_REF], originally developed for the high-dimensional context (see the related comments in the proof of Theorem IV.1 in this paper). In particular, the definitions [START_REF] Pipiras | Integral representations and of fractional Brownian motions[END_REF] and [START_REF]Exponents, symmetry groups and classification of operator frac-Brownian motions[END_REF] are applied in a mathematically sharp way. This establishes (ii).

Remark A.2. For a single j, we can show that the limit in distribution (31) of wavelet log-eigenvalues is scale-free. Consider the vector γ m as in (56) as well as the vector (60) In light of (48), it is clear that (60) converges to a Gaussian law that does not depend on j, as claimed.

To establish (60), recast entry-wise B(2 j ) =: b 1, 2 (2 j ) 1, 2=1,...,r (cf. (44)). By Proposition 3.1, (P5), in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF],

B(2 j ) = (2 j ) H 1 +H 2 +1 b 1, 2 (1)
1, 2=1,...,r .

(61)

Fig. 2 :

 2 Fig. 2: Repulsion bias in eigenvalue estimation. Top: Logarithms of estimated eigenvalues log 2 λm (2 j ) (left, blue) and log 2 λ m (2 j ) (right, red), as functions of the logarithms of the scales log 2 2 j = j for a synthetic mfBm with M = 6, N = 2 12 , Σ = I, W = I and H 1 = . . . = H 6 = 0.6, and j 2 = 8. Superimposed straight lines indicate the theoretical slopes 2H m +1 (top, black). Bottom: Logarithms of estimated eigenvalues log 2 λm (2 j ) (left, blue) and log 2 λ m (2 j ) (right, red) with theoretical linear behaviors subtracted.

Fig. 4 :

 4 Fig. 4: Estimation performance. Logarithm of the spectral norms of the bias (left), covariance (middle) and MSE (right) for the proposed bias-corrected multivariate estimator ĤM,bc (red), compared to those for the classical multivariate estimator ĤM (blue) and for the univariate estimator ĤU (black), as functions of the logarithm of the sample size N .

Fig. 5

 5 thus compares correlation coefficients between pairs ( Ĥm , Ĥm ), 1 ≤ m < m ≤ M , as functions of the sample size N , for the three different estimators ĤM,bc , ĤM and ĤU . The correlation structures are compared only for Config1 and Config3, in upper and lower triangles respectively. The conclusions drawn about the other configurations are similar. Fig. 5 clearly indicates that most pairs ( ĤU m , ĤU m ), 1 ≤ m < m ≤ M , remain significantly correlated, even for large sample sizes. By contrast, for the multivariate estimators, the pairs ( ĤM,bc m , ĤM,bc m ) and ( ĤM m , ĤM m ), 1 ≤ m < m ≤ M are quasi-decorrelated, for all sample sizes. These observations partly explain why the overall contribution of the variance to the MSE was larger for the univariate estimators when

Fig. 5 :

 5 Fig. 5: Correlation structures of ĤM,bc , ĤM and ĤU . Correlation coefficients between entries Ĥm and Ĥm , 1 ≤ m = m ≤ M , for ĤM,bc (red), ĤM (blue) and ĤU (black), as functions of sample sizes N , for Config1 (upper right triangle) and Config3 (lower left triangle) .

Fig. 6 :j 2 j

 62 Fig. 6: Approximation of the variance. Relative differences between empirical variances of the (left, black) univariate estimator ĤU m , (middle, blue) multivariate estimator ĤM m and (right, red) bias-corrected multivariate estimator ĤM,bc m and the approximate variance V N (Eq. (38)) with superimposed components m = 1, . . . , 6 as functions of the sample size N .

fww 2 ξ

 2 N,m (B) -f N,m (B(2 j )) = 1≤i1≤i2≤M ∂ ∂b i1,i2 f N,m ( B) • π i1,i2 B -B(2 j )(52)holds for some matrix B ∈ S >0 (M, R) lying in a segment connecting B and B(2 j ) across S >0 (M, R). Define the eventA N = {ω : B (τ )a (2 j ) ∈ O δ0 , τ = 1, . . . , 2 j 0 2 -j }. By (49), P(A N ) → 1, N → ∞. Hence, by (52), with probability going to 1, for large enough N , for τ = 1, . . . , 2 j 0 2 -j and m = 1, . . . , M , expansion (52) holds with B = B (τ ) a (2 j ), B = B(τ) a (2 j ) ∈ S >0 (M, R).In addition, can express, for anyi 1 , i 2 ∈ {1, . . . , M }, ∂ ∂b i1,i2f N,m ( B(τ) τ) (a2 j ) := Wa diag(H) B(τ) a (2 j )a diag(H) W * . Expressions (52) and (53) imply that, as N → ∞,ln 2 n a,j 0 2 f N,m ( B (τ ) a (2 j )) -f N,m (B a (2 j )) i1 (N ), u m (N ) a Hi 1 -Hm • w i2 (N ), u m (N ) a Hi 2 -Hm π i1,i2 B (τ ) a (2 j ) -B a (2 j ) i1 (N ), u m (N ) a Hi 1 -Hm • w i2 (N ), u m (N ) a Hi 2 -Hm π i1,i2 B (τ ) a (2 j ) -B(2 j ) .(54) However, by Proposition B.1 in[START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], under condition[START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF] we can suppose that there exists a deterministic unit vector u m such thatu m (N ) P → u m .(55)So, letγ m = {γ m,i } i=1,...,r := p-lim N →∞ W * u m (N ) ∈ R r .(56)Also, by Lemma B.8 in[START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], we can assume the existence of the limiting vectorx * ,m ={x * ,m,i } i∈I+ := p-lim N →∞ w i , u m (N ) • a(N ) Hi-Hm i∈I+ ∈ R r3 .(57) Therefore, (54) has the same limit in distribution as √ n a,j 0

m

  := B -1 33 (1)B 32 (1)π r1+1,...,r1+r2 (γ m ) ∈ R r3 . m (a(N )2 j ) -log 2 λ (τ ) m (a(N )2 j ) • m,i2 • π i1,i2 B (τ )a (1) -B(1) .

:

  Multivariate scale-free dynamics in EEG time for preictal and interictal states. Logarithms of the estimated diagonal entries log 2 S m,m (2 j ) (left, black) , estimated eigenvalues log 2 λm (2 j ) blue) and corrected eigenvalues log 2 m j red), the wavelet for m = 1, . . , for interictal (top) and preictal Equivalent conclusions were drawn for all other patients.

	Fig. 8: Differences in multivariate selfsimilarity between
	preictal and interictal states. Logarithms of the sorted p-
		for Wilcoxon rank-sum tests between preictal and
	interictal distributions for univariate ĤU m multivariate estimates ĤM m (blue) and bias
	ate estimates	(red), vs. Benjamini-Hochberg
	hypothesis testing thresholds dashed
	estimator	ĤM,bc	shows stronger in selfsimilarity
	parameters between preictal and interictal states, compared to
	either the estimates	ĤU	or the multivariate estimates
	ĤM	. This illustrates the benefits of taking into account
			dynamics in time series and
	correcting the bias multivariate estimation
	epilepsy seizure prediction.
	These analyses reported for single patient only (Patient
	22).

SUPPLEMENTARY MATERIAL

Fix

∈ N. Throughout the appendices, S( , R) and S ≥0 ( , R) denote, respectively, the space of symmetric and the cone of symmetric positive semidefinite × matrices. The notation vec S denotes the projection on the free entries of a symmetric matrix. Also, π i1,i2 and π i1,i2,...,i , ∈ N, denote the projection of a matrix or vector, respectively, on its entries i 1 , i 2 and i 1 , i 2 , . . . , i . We express any vector v entry-wise as {v } . Recall that n a,j = N a(N )2 j , j = j 0 1 , . . . , j 0 2 .

A. Proof of consistency

Proof. We only show (i). Fix j ∈ {j 0 1 , . . . , j 0 2 }, τ ∈ {1, . . . , 2 j 0 2 -j } and m ∈ {1, . . . , M }. Suppose, for the moment, that the limit [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF] holds. Recall that p-lim N →∞ denotes limit in probability. Writing a(N ) = a(N )2 j ,

This establishes the scaling relationship in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF].

We now need to show the convergence in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]. The argument follows from a simple adaptation of the technique for establishing Theorem 3.1 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], which originally pertains to a high-dimensional context.

In fact, in that theorem the measurements are given by a high-dimensional signal-plus-noise stochastic process of the form Y (t) = P(N )X(t) + Z(t), where, for p ≥ M , Y (t), Z(t) ∈ R p , P(N ) ∈ M(p, M, R) and X(t) ∈ R M . X(t) is the latent fractional stochastic process (the "signal"), Z(t) is "noise" component, P(N ) is a full rank of coordinates. So, first assume the p is fixed and set M . Now suppose X(t) an under the conditions of Theorem IV.1. Next set the matrix P(N ) to the M × M identity matrix. In addition, suppose the noise term Z(t) is identically zero. In order to establish [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], it now suffices to follow the argument of the proof of Theorem 3.1 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]. The convergence [START_REF]Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] can be established by an analogous argument. This shows (i).

Remark A.1. The framework for proving Theorem IV.1 provides a great deal of information about the asymptotic rescaled eigenvalues ξ m (2 j ), m ∈ {1, . . . , M }, as in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]. More specifically, for some function ϕ and some vector w, we can write p-lim

In order to describe ϕ( w), we need a few definitions and constructs. So, for the fixed m, define the sets of indices I -:= { : H < H m } and

Note that I -and I + are possibly empty. Also, let

be their respective cardinalities. For every τ = 1, . . . , 2 j 0 2 -j , define the auxiliary wavelet random matrix

) where a ≡ a(N ). Also, let B(2 j ) = EB a (2 j ), which neither depends on τ nor on a. Recast

In (44), B i and B i denote blocks of size r i ×r , and the scale 2 j is omitted wherever doing so is unambiguous. Similarly, define

where, for i = 1, 2, 3, W i ∈ M(p, r i , R), and each column vector is denoted by w , = 1, . . . , r.

Then, for ϕ as in (46), the proof of Proposition 5.1 and Lemma B.1 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF] imply that there exists a possibly random vector

based on which we can express the limit (40) as ϕ( w) (n.b.: ξ m (2 j ) is deterministic). In particular, relations (46) and (47) imply that ϕ(w) -and, hence, ξ m (2 j ) -only depends on the matrices B 22 , B 23 , B 33 , W 2 and W 3 . Consequently, the asymptotic rescaled function in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF] only depends on parameters associated with the index sets I 0 and I + in (41).

B. Proof of asymptotic normality

In the proof of Theorem IV.2, we write a ≡ a(N ) whenever notationally convenient. Also, we make use of the following notation. Let {X N } N ∈N , {Y N } N ∈N be two sequences of random variables. We write X N d ∼ Y N , N → ∞, when the two sequences have the same limit in distribution, which is assumed to be well defined.

Proof. First, we show (i). Fix j ∈ {j 0 1 , . . . , j 0 2 }. The proof of (31) is by means of an adaptation of the proof of Theorem 3.2 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]. In fact, let B a (2 j ) and B(2 j ) be as in (44). Under conditions (OFBM1-3), as a consequence of Theorem 3.1 in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF],

So, fix m ∈ {1, . . . , M }. Recall that, for any matrix M ∈ S(M, R), the differential of a simple eigenvalue λ m (M) exists in a connected vicinity of M and is given by

Bearing in mind (55) and (59), by the proofs of Lemma B.8 and Proposition 5.1 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], and by (61), we can recast (57) in the form

) Hm (2 j ) -diag(Hr 1 +r 2 +1,...,Hr ) m .

(62)

Rewrite, entry-wise,

By an adaptation of the proof of Proposition 3.1, (P6), in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF] we obtain

So, by (54), ( 61), ( 62) and (64), as N → ∞,

a (1) -B(1) . In view of the scaling relation in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], relation (60) holds, as claimed.

C. Approximation of eigenvalue covariance

Fix j ∈ {j 0 1 , . . . , j 0 2 }. In this appendix, we use the bivariate context to provide some justification for the very useful approximations

as well as

So, consider the case M = 2 and 0 < H 1 < H 2 < 1. Since the arguments for showing all relations (65)-(68) are similar, calculations are detailed for the bias-corrected eigenvalues λ(τ) m only. Let B a (2 j ) be as in (63) and further write entry-wise

In view of (60), it can be shown that we can write the expansions

and

(see relations (4.33) and (4.35) in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]). So, for the sake of computation, we normalize b 11 (1) = 1 = b 22 [START_REF] Frisch | Turbulence: the Legacy of A.N. Kolmogorov[END_REF]. Also, for notational simplicity, we write

where π 1 and π 2 denote the projection onto the -coordinate. For ofBm, the celebrated property of approximate decorrelation in the wavelet domain is given by Proposition 3.2 in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]. It implies that, for k = k ,

In view of (72), of the stationarity of fixed-scale wavelet coefficients (see [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], Proposition

(75) Namely, log 2 λ(τ) 1 (a(N )2 j ) and log 2 λ(τ) 2 (a(N )2 j ) are approximately decorrelated over large scales for large values of n a,j 0 2 . Alternatively, suppose condition (69) is dropped, namely, assume b 2 12 (1) = 0. Then, after appropriately eliminating null denominators, the expansions (70) and (71) and the approximation (73) immediately show that the log-eigenvalues are approximately decorrelated.

Next, we develop approximations for the variances. Similarly to (73) and ( 74


Marie Vayer

Robert Josef Mayer

Joseph Moran

David Leboeuf

Leveraging the Hydroarylation of α-(Trifluoromethyl)styrenes to Access Trifluoromethylated All-Carbon Quaternary Centers

Keywords: -(trifluoromethyl)styrenes, hydroarylation, hexafluoroisopropanol, Brønsted acid, all-carbon quaternary centers

published or not. The documents may come

In the last two decades, the trifluoromethyl group has emerged as a key motif in drug discovery, crop science and materials science due its ability to reshape the properties of a parent molecule, notably its permeability, lipophilicity, and metabolic stability. [START_REF] Müller | Fluorine in Pharmaceuticals: Looking Beyond Intuition[END_REF] Hence, several elegant methods have been designed to introduce this group in a straightforward fashion. [START_REF] Nie | Asymmetric Construction of Stereogenic Carbon Centers Featuring a Trifluoromethyl Group from Prochiral Trifluoromethylated Substrates[END_REF] Most of them rely on the installation of the trifluoromethyl group via different strategies such as nucleophilic, electrophilic, or radical trifluoromethylation. [START_REF] Nie | Asymmetric Construction of Stereogenic Carbon Centers Featuring a Trifluoromethyl Group from Prochiral Trifluoromethylated Substrates[END_REF] Recently, the direct functionalization of compounds in which the trifluoromethyl group is already pre-installed has become increasingly popular as a means of constructing trifluoromethylated all-carbon quaternary centers. An advantage of this approach is to rapidly build a large library of trifluoromethylated compounds from a single precursor by using different sets of carbon nucleophiles. Yet, to date, only a handful of transformations such as the 1,4-addition to activated alkenes, [START_REF] Kawai | Organocatalytic Asymmetric Synthesis of Trifluoromethylsubstituted Diarylpyrrolines: Enantioselective Conjugate Cyanation of β-Aryl-β-trifluoromethyl-disubstituted Enones[END_REF] the 1,6-addition to para-quinone methides, 4 the direct nucleophilic substitution of preactivated -trifluoromethyl alcohols, 5 and coupling reaction with allenes have effectively used this strategy (Scheme 1A). 6 An important limitation is that those methods often require tailored substrates that limit the scope of products that are attainable. 7 In light of this, it was surprising to find that the hydroarylation of unactivated -(trifluoromethyl)styrenes 8 was unprecedented, which might be explained by the relative robustness of the alkene functional group when compared to the trifluoromethyl group. Indeed, in a seminal study by the group of Ichikawa on Lewis acid-promoted arylation of -(trifluoromethyl)alkenes, only 3,3difluoroallylated arenes were obtained, 9 resulting from the activation of the CF3 group (Scheme 1B). Similar trends have also been observed under transition metal catalysis [START_REF] Miura | Synthesis of gem-Difluoroalkenes via β-Fluoride Elimination of Organorhodium(I)[END_REF] and photocatalysis, [START_REF] Wiles | Metal-Free Defluorinative Arylation of Trifluoromethyl Alkenes via Photoredox Catalysis[END_REF] as well as with organolithium reagents. [START_REF] Bégué | A Concise Synthesis of Functionalised gem-Difluoroalkenes, via the Addition of Organolithium Reagents to α-Trifluoromethylstyrene[END_REF] A hydroarylation process could facilitate the preparation of compounds incorporating trifluoromethylated all-carbon quaternary centers, offering at the same time more flexibility in terms of reaction scope since -(trifluoromethyl)styrenes can be rapidly prepared in one step by Suzuki-Miyaura cross-coupling or Wittig reaction from either commercially available 2-bromo-3,3,3trifluoro-1-propene or 2,2,2-trifluoroacetophenones (see Supporting Information for details). Such a reaction is, however, challenging. So far, there is only a single example of a hydrofunctionalization (hydroformylation) of -(trifluoromethyl)styrenes with the introduction of a functional group at the -position of the trifluoromethyl group. 7b Scheme 1. Strategies towards the preparation of trifluoromethylated all-carbon quaternary centers While the hydroarylation of unactivated styrenes has been extensively investigated, [START_REF] Rueping | A Review of New Developments in the Friedel-Crafts Alkylation[END_REF] the hydroarylation of highly electronically deactivated styrenes remains underdeveloped. However, we have recently demonstrated that the acidity of self-assembled H-bond clusters of hexafluoroisopropanol (HFIP) [START_REF]For selected reviews on HFIP[END_REF] could be amplified through coordination to a Lewis acid (Ca(NTf2)2/nBu4NPF6) or a Brønsted acid (TfOH) to trigger cationic reactions using highly deactivated substrates. 15,[START_REF] Vuković | For a related example on the Friedel-Crafts alkylation of highly electronically deactivated benzyl alcohols[END_REF] Inspired by those results, our plan was to use such a combination to enable the first catalytic intermolecular hydroarylation of -(trifluoromethyl)styrenes with readily available arenes as nucleophilic partners, thereby rapidly assembling trifluoromethylated diarylalkanes via the formation of a -(trifluoromethyl) carbenium. [START_REF] Fernandes | CF3-Substituted Carbocations: Underexploited Intermediates with Great Potential in Modern Synthetic Chemistry[END_REF] We began our optimization studies with 1,1,1trifluoromethylstyrene 1a as a model substrate in the presence of TfOH as a catalyst (10 mol%) in HFIP (0.2 M) as a solvent (Table 1). By using anisole 2a as an arene nucleophile, the target product 3 was obtained in 84% yield in the form of a separable 92:8 mixture of para/ortho regioisomers (entry 1). Reducing the temperature to 80 °C, the catalytic loading to 5 mol% or the number of equivalents of anisole resulted in decreased yields (entries 2-4). Increasing the concentration to 0.4 M proved to be detrimental to the reactivity (entry 5). Other Brønsted and Lewis acids such as HNTf2, [START_REF] Zhao | Triflimide (HNTf2) in Organic Synthesis[END_REF][START_REF]The reaction was also tested with other Brønsted acids such as HCl, p-TsOH, MsOH and chiral phosphoric acids. However, no reaction was observed, highlighting the critical role of TfOH as a catalyst[END_REF] Ca(NTf2)2/nBu4NPF6, 20 B(C6F5)3, 21 Re2O7 [START_REF] Volchkov | Recent Developments of Direct Rhenium-Catalyzed [1,3]-Transpositions of Allylic Alcohols and their Silyl Ethers[END_REF] and Bi(OTf)3 [START_REF] Ondet | Cyclisations Catalysed by Bismuth(III) Triflate[END_REF] were tested but did not react or generate 3 in a synthetically useful yield (entries 6-10). On the other hand, the combination of Bi(OTf)3/nBu4NPF6 afforded 3 in 71% yield (entry 11).

Table 1. Optimization of reaction conditions

Entry

Variation from standard conditions [a] Yield [b] (%) (p/o ratio) 1,2-DCE instead of HFIP 7 (82:18) [d] 14

MeNO2 instead of HFIP - [e] 15 toluene instead of HFIP traces 16 without TfOH n.r.

[a] Standard reaction conditions: 1a (0.2 mmol), anisole (5.0 equiv.), TfOH (10 mol%) in HFIP (0.2 M), 110 °C, 24 h (in a sealed tube). [b] Isolated yield. [c] 5.0 mmol scale. [d] NMR yield using triethylsilane as an external standard. [e] A mixture of several products was observed.

Other common organic solvents were unable to provide 3 (entries 12-15), highlighting the crucial role of HFIP in this reaction. Finally, a control reaction in the absence of TfOH was conducted, and no product was observed (entry 16).

From a practical point of view, upscaling the reaction (5.0 mmol) allowed a gram-scale synthesis of the trifluoromethylated compound 3 in 79% yield (3.96 mmol, 1.11 g) (entry 1).

With those reaction conditions in hand, we proceeded to investigate the scope of this reaction, initially focusing on the reactivity of various electron-rich arenes with 1a (Scheme 2). In most cases, the reaction was achieved in high Scheme 2. Scope of arenes 2a-2u with 1,1,1 trifluoromethylstyrene 1a

[a] Isolated yield of the major regioisomer. [b] NMR yield and ratio p/o using triethylsilane as an external standard. [c] Reaction temperature 80 °C.

yields (up to 96%), along with an excellent regioselectivity. Toluene (4), cumene (5), naphthalene (6), 24 biphenyl (7), xylene (8-9), tetrahydronaphthalene (10), phenol (11-15) and sulfide (15) derivatives were among the arene nucleophiles compatible with our standard conditions. Moreover, the combination of a donor substituent (OMe) with a moderate electron-withdrawing group (Cl) yielded the corresponding products in high yields (16-18, 70-95%).

In several cases, regioisomers could be separated by flash column chromatography. Importantly, less nucleophilic arenes such as benzene (20, 67%) also proved to be competent nucleophiles for this transformation. In turn, an electron-poor (bromobenzene) or a sterically hindered arene (mesitylene) reacted poorly under the reaction conditions, leading to a slow decay of the substrate 1a over time. The same observation was made with when indole was used as nucleophile.

Next, we examined the scope of trifluoromethylated styrenes using anisole 2a, 1,2-chloroanisole 2o, and thiophene 2v as nucleophiles (Scheme 3). The reaction was compatible with electron-donating substituents at the para Scheme 3. Scope of 1,1,1-trifluoromethylated styrenes 1a-1k

[a] Isolated yield. [b] NMR yield using triethylsilane as an external standard. [b] Reaction temperature 80 °C. [c] Reaction time 48 h in the presence of 20 mol% of TfOH.

and meta position, delivering products 21-25 in yields ranging from 33 to 97%. In some cases (22 and 23), products bearing two electron-rich arenes were prone to retro-hydroarylation upon purification by flash column chromatography, leading to diminished yields. Other electron-rich (hetero)arenes such as biphenyl (26, 73%), naphthalene (27, 83%) and benzothiophene (28, 49%) were also amenable to the reaction, as well as moderate electron-poor arenes such as aryl bromides (29-31, 12-82%). However, in the case of substrates incorporating an arene bearing an electron-withdrawing group, such as an ester (32), the substrate remained intact as it is likely too electronically deactivated to react. We then turned our attention to non-terminal alkenes. In the case of indene 1j, trifluoromethylated indanes 33 and 34 were obtained in high yields (88% and 87%, respectively). 25 In turn, an acyclic trisubstituted alkene also provided the target product 35, albeit in a moderate yield (33%), which is the result of a competitive hydroarylation at the other position of the alkene (30% yield).

Of note, the presence of an ortho-OMe substituent posed complications, as the hydroarylation becomes reversible under the reaction conditions. For example, using 2o as arene nucleophile led to a 42% yield of the re-arranged product 22, arising from a retro-hydroarylation of the omethoxyphenyl moiety and its subsequent alkylation at the para position with the newly formed styrene 37 (Scheme 4).

Scheme 4. Retrohydroarylation with substrate 1l

Importantly, this transformation is not limited to a CF3 group but could be extended to the CF2H group, [START_REF] Sap | Late-Stage Difluoromethylation: Concepts, Developments and Perspective[END_REF] whose properties differ from the CF3 group. [START_REF] Zafrani | Difluoromethyl Bioisostere: Examining "Lipophilic Hydrogen Bond Donor" Concept[END_REF] In the case of substrate 1m, the reaction with anisole led to the product 38 in 87% (para/ortho 86:14). It is noteworthy that the difluoromethyl analog is far more reactive as the reaction was achieved in only 1 h at 60 °C.

To illustrate the synthetic utility of this transformation, several derivatizations were carried out to introduce functional groups of interest (Scheme 5). For instance, the oxidation of 1,4-benzenediol 14 was conducted in the presence of cerium ammonium nitrate (CAN) to deliver quinone 39 in a nearly quantitative yield. In addition, the regioselective formylation of phenol 11 was achieved in 72% yield (40) by using the Duff reaction. [START_REF] Lindoy | Mono-and Diformylation of 4-Substituted Phenols: A New Application of the Duff Reaction[END_REF] As mentioned above, the hydroarylation reaction was not compatible with substrates that are deactivated by the presence of a second electron-poor functional group; however, this issue may be circumvented by a palladium-catalyzed cross-coupling reaction between bromobenzene 29 and zinc(II) cyanide (41, 64% yield). In another application, the pmethoxyphenyl group of 3 can be used as a masked carboxylic acid via oxidation (42, 46%), [START_REF] Carlsen | A Greatly Improved Procedure for Ruthenium Tetroxide Catayzed Oxidations of Organic Compounds[END_REF] thereby mimicking the hydrocarboxylation of 1,1,1trifluoromethylstyrene 1a.

Scheme 5. Derivatizations of trifluoromethylated compounds

To gain insights into the mechanism of the reaction, DFT computations were performed on the reaction of 1a with benzene at the SMD(HFIP)/MN15/def2-TZVP level of theory (Figure 1). 30 Protonation of 1a by HOTf can either yield the carbenium ion 1a-H + or, via significantly less favorable HF elimination, the difluorocarbenium ion A. The large energy difference computed for the formation of 1a-H + and A is in agreement with our experimental findings that only products originating from + experimentally obtained and not those resulting from HF elimination via cation A. Though formation of cation A 1a-H is unfavorable our conditions, stoichiometric addition of a Lewis acid that binds fluoride can provide the driving force for the selective formation of A reported by Ichikawa. 9a Once 1a-H + can react multiple available nucleophiles in solution. The reactions with triflate or HFIP to yield the trapping products 1a-OTf and 1a-HFIP, respectively, are computed to be only slightly exergonic and are thus expected to be reversible. 31 Nucleophilic addition of benzene to the carbenium ion 1a-H + proceeds via a moderate activation barrier relative to the carbenium ion (+48.9 kJ mol -1 ) to yield the Wheland intermediate. Deprotonation of this intermediate by any available base regenerates the acid catalyst and yields the final product 20, thus closing the catalytic cycle.

Finally, we compared the cation 1a-H + with the cumyl cation C. Starting from HOTf and either 1a or -methylstyrene, the Gibbs free energies for the formation of the carbenium ions are +72.6 or -3.2 kJ mol -1 , respectively, suggesting that the formation of 1a-H + is significantly more endergonic. However, once formed 1a-H + is a significantly stronger electrophile then C: In the reaction with benzene, 1a-H + reacts via a 32.6 kJ mol -1 lower activation barrier then cation C (for the full Gibbs energy profile starting from methylstyrene, see Supporting Information, Figure S2). The different stability of 1a-H + and C is furthermore illustrated by gas-phase methyl anion affinities (GMA, calculated at the B3LYP/6-311++G(3df,2pd)//B3LYP/6-31G (d,p) level). Methyl anion affinities are an established tool to compare cation stabilities and have previously been computed for a wide range of benzhydrylium ions that are standard tools in carbocation chemistry. [START_REF] Mayr | Scales of Lewis Basicities toward C-Centered Lewis Acids (Carbocations)[END_REF] While the cumyl cation is found to have a similar GMA compared to the parent benzhydrylium Figure 1. (A) Gibbs energy profile calculated at the SMD(HFIP)/MN15/def2-TZVP level of theory and (B) comparison of the activation barriers for benzene addition of 1a-H + and the cumyl cation C at the same level. (C) Gas-phase methyl anion affinities (GMA) calculated at the B3LYP/6-311++G(3df,2pd)//B3LYP/6-31G (d,p) level. a When using one explicit HFIP molecule on both TfOH and TfO -, an insignificantly different value of +73.6 kJ mol -1 was obtained for the formation of 1a-H + starting from 1a. ion Ph2CH + (GMA = -827.3 kJ mol -1 ), the value obtained for cation 1a-H + is higher than for any of the benzhydrylium ions previously characterized (see Supporting Information, Figure S3 for a comparison). [START_REF] Mayr | Scales of Lewis Basicities toward C-Centered Lewis Acids (Carbocations)[END_REF] This furthermore illustrates the unusually high instability and thus reactivity of cation 1a-H + compared to more classical carbenium ions and the abilities of HFIP to enable its formation. One potential role of HFIP in mediating this transformation might again be in reversibly trapping the carbenium ion 1a-H + e.g. in the form of 1a-HFIP, thus preventing potential oligomerization of 1a observed in its absence. 30 In conclusion, we herein presented a new hydroarylation reaction that directly forms trifluoromethylated all-carbon quaternary centers. This method relies on the unique combination of TfOH and HFIP to promote the target reaction, whose mechanism was rationalized by DFT computations. This reaction is facile to operate, scalable and provides most of the products in high yields and with good regiocontrol. The compounds obtained can serve as a platform to access more complex molecules. Therefore, we anticipate that this new route to trifluoromethylated compounds will open new opportunities in drug discovery.

ASSOCIATED CONTENT

Supporting Information.

The Supporting Information is available free of charge via the Internet at http://pubs.acs.org. Experimental procedures, characterization data, copies of NMR spectra, procedure of quantum-chemical calculations and xyz coordinates of all optimized structures (PDF)

1 none 84 (92:8), 79 [c]

	2	80 °C instead of 110 °C	73 (88:12)
	3	5 mol% of TfOH instead of 10 mol%	55 (86:14)
	4	3 equiv. of 2a instead of 5 equiv.	60 (92:8)
	5	0.4 M instead of 0.2 M	64 (92:8)
	6	HNTf2 instead of TfOH	16 (87:13) [d]
	7	Ca(NTf2)2/nBu4NPF6 instead of TfOH	5 (87:13) [d]
	8	B(C6F5)3 instead of TfOH	n.r.
	9	Re2O7 instead of TfOH	traces
	10	Bi(OTf)3 instead of TfOH	33 (87:13)
	11	Bi(OTf)3/nBu4NPF6 instead of TfOH	71 (87:13)
	12	TFE instead of HFIP	5 (83:17) [d]
	13

ACKNOWLEDGMENT

This work was supported by the Interdisciplinary Thematic Institute ITI-CSC via the IdEx Unistra (ANR-10-IDEX-0002) within the program Investissement d'Avenir. R.J.M. thanks the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for a fellowship (MA 9687/1-1). Computations were performed on the High Performance Computing Center of the University of Strasbourg. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data. D.L. and J.M thank the CNRS. The authors thank Wahnyalo Kazone for HRMS analysis.

AUTHOR INFORMATION

Corresponding Author dleboeuf@unistra.fr

Author Contributions

The manuscript was written through contributions of all authors.

Notes

The authors declare no competing financial interest.