Davide Faranda

Bérengère Dubrulle

François Daviaud

Flavio Maria

Emanuele Pons

Franc ¸ois Daviaud

Maria Emanuele Flavio

Pons

Probing turbulence intermittency via autoregressive moving-average models

Keywords: number(s): 47, 27, E

Introduction. One of the few exact results known for isotropic, homogeneous, and mirror-symmetric turbulence is the 4/5 law derived by Kolmogorov. It links the longitudinal velocity increments δu = u(x + )u(x) to the mean rate of energy dissipation via

δu 3 = -4 5 , (1)
where angular brackets denote averaging. This exact relation was generalized by Kolmogorov [1] as a scaling law δu ≡ ( ) 1/3 , where ≡ means has the same statistical properties. Should be a nonstochastic constant, the scaling law would imply self-similar behavior for the structure functions of order p, S p ( ) = δu p , which would scale like F p ( ) ∼ p/3 p/3 .

(

) 2
For p = 3, we recover the 4/5 law. For p = 2, this equation predicts a second-order structure function that varies like 2/3 . By a Fourier transform, this is equivalent to a one-dimensional energy spectrum scaling with wave number k as E(k) ∼ k -5/3 , also known as the Kolmogorov spectrum [2,3]. Both the 4/5 law and the Kolmogorov spectrum have been measured and checked in many natural and laboratory isotropic turbulent flows [START_REF] Frisch | Turbulence[END_REF]. More generally, Eq. ( 2) predicts a linear law for the exponent of the structure functions ζ (p) = d ln F p ( )/d ln = p/3. However, as pointed out by Landau and recognized by Kolmogorov [1], there is no reason to assume that is a constant over space and/or time, so it should rather be viewed as a stochastic process, which depends upon the scale at which it is measured ≡ ( ). In such a case, the correct scaling of the structure function is rather

F p ( ) ∼ ( ) p/3 p/3 . ( 3
)
This modified law predicts correction to the linear law ζ (p) = p/3, which is connected to the intermittent nature of the dissipation. For example, a log-normal model for the dissipation (a suggestion by Landau and Obukhov) implies * davide.faranda@cea.fr quadratic corrections for the ζ (p). Other models have been suggested and lead to different corrections [START_REF] She | [END_REF][6][7]. Intermittency corrections up to p = 4 have been measured in a variety of experimental and numerical flows and appear to be robustly consistent from one experiment to another (see, e.g., the review of [8]). Corrections for larger values of p are subject to resolution and statistical convergence issues: The larger the scaling exponent, the larger the statistical sampling must be in order to capture the rare events. There is presently no general consensus about the behavior of intermittency corrections at large order. This hinders progress in the understanding of the statistical properties of the energy dissipation. In this Rapid Communication we suggest an approach to probing intermittency corrections based on the autoregressive movingaverage (ARMA) modeling of turbulent time series. We introduce an index ϒ that measures the distance from a Kolmogorov-Obukhov model in the ARMA space. Applying our analysis to velocity measurements in a von Kármán swirling flow, we show that this index is proportional to the traditional intermittency correction computed from the structure function and provides the same information, using shorter time series. Intermittency parameters. In most laboratory turbulent flows, data sets are time series of values of a physical observable at a fixed point or obtained by tracking Lagrangian particles. This motivated the shift of paradigm from space velocity increments to time velocity increments defined as δu τ = u(t + τ )u(t) and motivated measurements of time structure functions G p (τ ) = (δu τ ) p and its local exponent χ p = d ln G p (τ )/d ln τ . In situations where measurements are made on the background of a strong mean velocity U , scale velocity increments and time velocity increments can be directly related through the Taylor hypothesis = Uτ . When fluctuations are of the same order as the mean flow, however, the Taylor hypothesis fails. A suggestion has been made in [9] to then resort to a local Taylor hypothesis, in which = dtu(t), where u is the local rms velocity. This is equivalent to considering a scale such that ∼ τ δu τ and may be seen as equivalent to modifying the space-refined Kolmogorov hypothesis into a time hypothesis δu τ ≡ ( τ ) 1/2 , which leads to

G p (τ ) ∼ p/2 τ p/2 . ( 4
)
Such scaling is equivalent to the scaling obtained using the Lagrangian structure function. In any case, we may define the intermittency as the deviation of the local exponents ζ * p = ζ p (space increments) or ζ * p = χ p (time increments) with respect to a linear behavior and may be quantified to first order by the parameter

μ = ζ * 2 -2 3 ζ * 3 . ( 5
)
This factor is proportional to the logarithm of the β parameter of the log-Poisson model [START_REF] She | [END_REF]6], or to the μ parameter of the lognormal model [1]. It is also valid when the scaling exponents have been computed using extended self-similarity (ESS) [10], which is especially interesting in situations where turbulence is inhomogeneous and when the Taylor hypothesis does not hold. In the following, we compare this intermittency index with another one, built in a purely statistical framework. Indeed, Thomson [11] showed that, in the Lagrangian framework, the time-refined Kolmogorov hypothesis is in fact equivalent to a stochastic description in terms of an Ornstein-Uhlenbeck process with suitable drift and noise term

du = - u T dt + C 0 dW, ( 6
)
where T is a decorrelation time scale, C 0 is a universal constant, and is the mean dissipation. Indeed, taking into account the definition of the particle position x, dx = udt, we get a scaling of the time averages of velocity and position as

u 2 (t) ∼ t, x 2 (t) ∼ t 3 . ( 7
)
The second property is the Richardson law. Then, defining δu = [u 2 (t)] 1/2 and = [x 2 (t)] 1/2 , we get from Eq. ( 7) δu ∼ 1/3 , which leads to the space-refined Kolmogorov hypothesis. The discrete-time version of Eq. ( 6) can be written as

u t = φu t-1 + ψ t , ( 8
)
where t is a discrete-time label, dW are the increments of a Brownian motion, φ = (1 -t T ), and ψ t are independent variables, normally distributed. Equation ( 8) is the expression of an autoregressive process of order one, denoted by AR (1). Such a model is described by a single decorrelation time. So it cannot describe real flows intermittency, which involves a whole range of time scales corresponding to the turnover times of the turbulent eddies with memory effects. To capture these effects, it is mandatory to consider a projection of the velocity data on higher-order ARMA(p,q) models. This enables a quantification of intermittency effects as a distance with respect to the insufficient AR(1) model in this space.

Intermittency as a distance in ARMA space. A stationary time series X t is said to follow an ARMA(p,q) process if it satisfies the discrete equation

X t = p i=1 φ i X t-i + ε t + q j =1 θ j ε t-j , (9)
with ε t ∼ N W (0,σ 2 ), where N W stands for white noise and the polynomials φ(z

) = 1 -φ 1 z t-1 -• • • -φ p z t-p and θ (z) = 1 -θ 1 z t-1 -• • • -θ q z t-q , with z ∈ C, have no common
factors. Notice that the noise term ε t will be assumed to be a white noise, which is a general condition [START_REF] Brockwell | Time Series: Theory and Methods[END_REF]. We ensure unicity by applying the Box-Jenkis procedure [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF]: We choose the lowest p and q such that the residuals of the series filtered by the process ARMA(p,q) are not correlated. To define a suitable distance in the space of ARMA(p,q) models, we introduce the Bayesian information criterion C BI , measuring the relative quality of a statistical model, as

C BI = -2 ln L(n, σ 2 ,p,q) + k[ln(n) + ln(2π )], (10)
where L(n, σ 2 ,p,q) is the likelihood function for the investigated model, k = p + q, and n is the length of the sample. The variance σ 2 is computed from the sample and is a series-specific quantity. The normalized distance between the fit ARMA(p + 1,q) and the Kolmogorov AR(1) model is then defined as the normalized difference between the C BI (n, σ 2 ,p + 1,q) and the AR(1) C BI (n, σ 2 ,1,0):

ϒ = 1 -exp{|C BI (p + 1,q) -C BI (1,0)|}/n. ( 11
)
The p + 1 serves to magnify ϒ near zero. Note that 0 ϒ 1: It goes to zero if the data set is well described by an AR(1) model and tends to one in the opposite case. In the case of velocity increments time series, it measures deviations from the Kolmogorov model. Application to turbulent data. We apply the index defined in Eq. ( 11) to velocity time series obtained in a von Kármán turbulent swirling flow. The experimental setup consists of two sets of blades mounted on two counterrotating coaxial impellers at the top and bottom of a cylindric vessel of diameter R = 0.1 m. The operating fluid is water and the rotation frequency of the impellers can reach F = 15 Hz, resulting in large Reynolds numbers (Re = 2πF R 2 ν -1 ∼ 10 6 ). A detailed description of the experiment can be found in [START_REF] Ravelet | [END_REF][15][16]. Two techniques are used to measure the fluid velocity on a grid: the particle interferometry velocimetry (PIV) and the laser Doppler velocimetry (LDV), mapped on a regular sampling time applying a sample-and-hold algorithm. The stereoscopic PIV measures the three components of the velocity field in a plane, while the LDV measurements provide the out-of-plane velocity component V φ in a plane. The PIV produces regularly sampled time series at intervals of 0.1 s over a sample size at most of order 10 4 and a spatial resolution of the order of 1 mm, i.e., 10-100 times larger than the dissipation scale. The LDV time series are sampled over a time scale of the order of 0.001 s, producing a sample size up to 10 6 data on a grid of spatial resolution of the order of 1 cm. Given these resolution constraints, we compute spatial (temporal) velocity increments for the PIV (LDV) data. The idea is to compute at each spatial grid location the classical intermittency index μ, compare it to ϒ, and see how they vary. All the analyses presented in this Rapid Communication are done using three components for the PIV and V φ for the LDV. Since the von Kármán flow is inhomogeneous and anisotropic with large fluctuations [16], we expect that the time and space velocity structure functions depend on the measurement points. This is illustrated in Fig. 1 for the second-and fourth-order spatial and time structure functions. For the spatial case, deviations from the Kolmogorov scaling (solid lines) are small for the spatial structure functions, near the symmetry plane Z = 0. This plane is the location of an intense shear layer and has traditionally been used to perform isotropic homogeneouslike measurements. Outside this plane, deviations from the Kolmogorov scaling are large. For the time case, one observes two distinct behaviors: Outside the shear layer, where a mean velocity is well defined, one observes close to Eulerian-Kolmogorov scaling at the smallest time increments τ p/3 ; in the shear layer, where no Taylor hypothesis holds, the scaling is closer to Lagrangian scaling τ p/2 . However, as already noted by [9] and shown in Fig. 2, the relative scaling exponents

ζ * p computed as G p (τ ) ∼ |δu τ | 3 ζ * p (ESS method
) are in most of the flow close to the universal scaling exponents found by [8], in a variety of homogeneous turbulent flows, even those with no obvious inertial range. Using these ESS scaling exponents to compute the μ index, we may then draw a map of the intermittency and compare it with ϒ. This is done in Fig. 3 for an LDV experiment at Re ∼ 10 5 . The spatial patterns look indeed similar. Moreover, the plot of ϒ as a function of μ [Fig. 3(c)] evidences a linear relation between them; the linear regression represented by the red line leads to a linear correlation coefficient r 0.69. This means that ϒ traces the same intermittency characteristics as the time structure functions. The comparison of ϒ with the intermittency index μ computed for spatial structure functions is also informative: Because of convergency issues, we have to use a data set p computed for the LDV experiments. Different lines correspond to different measure points. Red (gray) spots mark the scaling exponents reported in [8]. of about 10 5 -10 6 data points to converge the estimate of μ, while only 10 3 are needed to converge ϒ. To illustrate this, we use the longest data set available: 9000 velocity fields of a PIV experiment performed at Re 5 × 10 4 . At this value, the von Kármán flow experiences the equivalent of a phase transition [17], with time wandering of the shear layer in between Z = 0.3 and -0.3. This corresponds to a very large time intermittency and is detected by the ϒ index as shown in Fig. 4, under the shape of two patches at R 0, Z = 0.3 and R 0, Z = -0.3. This pattern is unique to the phase transition and is not present in other PIV experiments [18]. Besides, one observes a fairly symmetric structure, with maxima corresponding to the flow's four-cell structure structure functions, resulting in a lack of symmetry of the μ field (not shown). As a result, μ fluctuates over a decade around a value of about 0.05, while ϒ spans several orders of magnitude, as can be seen in Fig. 4(b). This shows that ϒ is a more sensitive tool to detect intermittency than μ.

Discussion. We have introduced an intermittency index ϒ that can be interpreted as a statistical distance between the best-fit linear ARMA model for a turbulent time series and the simplest possible process, i.e., AR(1). We have compared such an index with a classical intermittency index μ = ζ * 2 -2 3 ζ * 3 . In statistically converged LDV series, the two parameters are linearly related, with a regression coefficient R 0.69. In shorter PIV time series, ϒ catches important characteristics of the mean flow in situations where μ cannot be computed due to a lack of convergence. Therefore, the main advantage of this index is the applicability to cases in which no big data sets are available. Moreover, ARMA models contain the information on nonlocal interactions between large and small scales in the parameters φ i and θ i . In general, the AR(p) part of an ARMA(p,q) process is the contribution of the large scales and represents the persistence of the process. If the process is more persistent than AR(1) with |φ| < 1, higher values of p are required to explain all the correlation coming from the large scales. Analogous considerations hold for the MA(q) part. When there is a clear separation between large and small scales. On the other hand, when their magnitude is similar such separation is not ensured anymore [18]. Our results reinforce the hypothesis of [19] that intermittency propagates in direct interactions between large and small scales, rather than in cascades. Finally, our method can be used to validate models based on the stochastic differential equation [20] with respect to experimental data.

FIG. 1 .

 1 FIG. 1. (Color online) (a) Space structure function F p ( ) of order 2 in black and of order 4 in blue (gray) at two PIV grid points of coordinates R = -0.10, Z = -0.14 (circles) and R = 0.98, Z = 0.70 (crosses). (b) Time structure function G p (τ ) of order 2 in black and of order 4 in blue (gray) at two LDV grid points of coordinates R = 0, Z = 0.35 (circles) and R = 0.52, Z = 0.61 (crosses). Lines represent the Kolmogorov predictions: solid is Eulerian with 2/3 and τ 2/3 in black and 4/3 or τ 4/3 in blue (gray) and dotted is Lagrangian with τ in black and τ 2 in blue (gray).

FIG. 2 .

 2 FIG. 2. (Color online) Plot of ζ *p computed for the LDV experiments. Different lines correspond to different measure points. Red (gray) spots mark the scaling exponents reported in[8].

FIG. 3 .FIG. 4 .

 34 FIG. 3. (Color online) (a) Index ϒ vs (b) the intermittency index μ = ζ * 2 -2 3 ζ * 3 . White crosses show measurement points. (c) Scatter plot of ϒ vs μ. The red (gray) line shows a linear regression of the data.

Acknowledgments. Experimental data were provided with the help of E. Herbert, P. Ph. Cortet, C. Wiertel, and V. Padilla. D.F. acknowledges the support of a Centre National de la Recherche Scientifique postdoctoral grant.


Davide Faranda
email: davide.faranda@cea.fr

) Flavio

Maria Emanuele Pons

B Éreng Ère Dubrulle

Franc ¸ois Daviaud

Brice Saint-Michel

Éric Herbert

Pierre-Philippe Cortet

Modelling and analysis of turbulent datasets using Auto Regressive Moving Average processes

We introduce a novel way to extract information from turbulent datasets by applying an Auto Regressive Moving Average (ARMA) statistical analysis. Such analysis goes well beyond the analysis of the mean flow and of the fluctuations and links the behavior of the recorded time series to a discrete version of a stochastic differential equation which is able to describe the correlation structure in the dataset. We introduce a new index ϒ that measures the difference between the resulting analysis and the Obukhov model of turbulence, the simplest stochastic model reproducing both Richardson law and the Kolmogorov spectrum. We test the method on datasets measured in a von Kármán swirling flow experiment. We found that the ARMA analysis is well correlated with spatial structures of the flow, and can discriminate between two different flows with comparable mean velocities, obtained by changing the forcing. Moreover, we show that the ϒ is highest in regions where shear layer vortices are present, thereby establishing a link between deviations from the Kolmogorov model and coherent structures. These deviations are consistent with the ones observed by computing the Hurst exponents for the same time series. We show that some salient features of the analysis are preserved when considering global instead of local observables. Finally, we analyze flow configurations with multistability features where the ARMA technique is efficient in discriminating different stability branches of the system. C 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4896637]

I. INTRODUCTION

For a long time, experimentally testable predictions of turbulence properties have been influenced by available measurements. For example, hotwire velocity measurements have motivated statistical analysis of turbulent spectra or velocity increments computation, allowing the evaluation of Kolomogorov direct or refined similarity hypothesis. [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers[END_REF] More recently, new sophisticated instruments and acquisition techniques, such as the Particle Image Velocimetry (PIV) and the Laser Doppler velocimetry (LDV), have made possible to measure instantaneous velocity fields with resolution equivalent to that of Large Eddy Simulations. [START_REF] Kato | Laser doppler velocimeter[END_REF][START_REF] Adrian | Particle Image Velocimetry[END_REF][START_REF] Goldstein | Fluid Mechanics Measurements[END_REF] With these high quality datasets, it is now possible to reconstruct the large scale flow dynamics and compute global observables even in relatively complex geometries such as in the non-homogenous, non-isotropic von Kármán flow. [START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF][START_REF] Cortet | Normalized kinetic energy as a hydrodynamical global quantity for inhomogeneous anisotropic turbulence[END_REF] As more spatial and temporal scales are becoming accessible to measurements, it is important to extract all the possible information from a statistical analysis of the data as it may eventually lead to test new theoretical predictions. In that respect, the integrated information obtained by measuring spectral features may not be sufficient to distinguish the contributions of different scales. Moreover, filter response functions used in spectral analysis may introduce spurious effects on the fast scales hiding an intricate structure. [START_REF] Zhou | Interacting scales and energy transfer in isotropic turbulence[END_REF] Since present accessible measurements now give access to a large range of spatial scales, it therefore seems more promising to focus on turbulence properties in the physical space.

Alternative statistical description of turbulence in the physical space actually date back to Kolmogorov and Obukhov [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers[END_REF][START_REF] Obukhov | On the distribution of energy in the spectrum of turbulent flow[END_REF] and motivated formulation of stochastic models for the time evolution of turbulent observables. A now classical example is the Lagrangian stochastic model for the velocity of a passive tracer proposed by Thomson. 9 In that model, the inertial range of Lagrangian turbulent velocity is described through a Langevin equation, involving parameters that are determined via the so-called Well Mixed Condition (WMC). This model is in fact equivalent to an autoregressive process of order 1, usually denoted AR (1) or ARMA(1, 0) (see below) (ARMA = Auto Regressive Moving Average). Recent experiments however indicate that this simple model does not work for the velocity increments, which cannot be described by a simple standard Brownian motion as suggested by Obukhov. [START_REF] Frisch | Turbulence[END_REF] Indeed, non normal corrections originate from long correlations due to the intermittent character of turbulent flows. There are several models that suggest a more refined description, based, e.g, on the Rapid Distorsion Theory, [START_REF] Laval | Langevin models of turbulence: Renormalization group, distant interaction algorithms or rapid distortion theory?[END_REF] on the account of the two-point two-time Eulerian acceleration-acceleration correlation, [START_REF] Friedrich | Statistics of lagrangian velocities in turbulent flows[END_REF] on temporal memory kernel. [START_REF] Baule | Joint probability distributions for a class of non-Markovian processes[END_REF] These approaches lead to excellent approximations of the experimentally determined velocity pdf's, although an analytic solution for the model is still not available (for a review see Refs. 14 and 15). However, it is not clear whether these models directly correspond to the features really observed in turbulent experiments or, in other words, how far is a real experiment from the theoretical idealization.

To answer this question, as well as optimizing the information available from experimental measurements, it is mandatory to consider a more refined statistical analysis, able to account for temporal memory effects as well as velocity dependent diffusion coefficients. A good candidate is given by analysis in terms of ARMA(p,q) processes that have already been used to study problems ranging from geophysics to social science and finance. [START_REF] Franses | Non-linear Time Series Models in Empirical Finance[END_REF][START_REF] Guiot | ARMA techniques for modelling tree-ring response to climate and for reconstructing variations of paleoclimates[END_REF] This analysis aims to represent the statistical properties of a time series X t using a model in which the value at time t is a combination of the p previous observations of the series-the so-called auto-regressive part AR(p) -and q noise terms -the moving average part MA(q)-with p and q chosen to be the lowest order to describe the series (see below). We observe that ARMA(p, q) processes are also good candidates to describe turbulent experimental data, since high p orders correspond to high temporal memory and high q orders correspond to a complicated structure of the diffusion coefficients. In the present paper, we will apply the ARMA modeling technique to large datasets obtained in the (inhomogeneous, anisotropic) von Kármán flow to illustrate the potential of this method.

The von Kármán experiment, in which the flow is generated in between two counter-rotating coaxial impellers, is a simple way to obtain experimentally a large Reynolds number (Re ∼ 10 6 ) in a compact design. [START_REF] Marié | Experimental measurement of the scale-by-scale momentum transport budget in a turbulent shear flow[END_REF] In the equatorial shear layer, fluctuations are large and exhibit similar local properties as in large Reynolds number experimental facilities devoted to homogeneous turbulence. Away from the shear layer, one observes a decrease of the turbulence intensity. Overall, the flow is strongly turbulent, so that the instantaneous velocity fields, measured by means of a PIV system, strongly differ in a non-trivial manner from their time average. [START_REF] Cortet | Normalized kinetic energy as a hydrodynamical global quantity for inhomogeneous anisotropic turbulence[END_REF] Although significant advancements in understanding the physics of this system by statistical analysis [START_REF] Cortet | Experimental evidence of a phase transition in a closed turbulent flow[END_REF] and from statistical mechanics approaches [START_REF] Naso | Statistical mechanics of Beltrami flows in axisymmetric geometry: Theory reexamined[END_REF] have been made, several features of the flow remain unexplained and require further investigations. These include the nature of the phase transition recently discovered in the fully turbulent regime, [START_REF] Cortet | Experimental evidence of a phase transition in a closed turbulent flow[END_REF] the forcing dependent stability of steady states [START_REF] Saint-Michel | Evidence for forcing-dependent steady states in a turbulent swirling flow[END_REF] or the asymmetry of the torque probability distribution in different forcing conditions. [START_REF] Leprovost | A stochastic model of torques in von kármán swirling flow[END_REF][START_REF] Titon | The statistics of power injected in a closed turbulent flow: Constant torque forcing versus constant velocity forcing[END_REF] These features are based on both local measurements (such as velocity measurements using PIV or LDV techniques) or global measurements, such as total angular momentum, energy or torque applied to the rotating disks by the turbulence (drag friction). For any of these local and global measurements, we will define the ARMA(p, q) model which better represents the data, keeping in mind that the simplest model explaining the Kolmogorov turbulent spectra is the ARMA(1,0) model (see below). We then try to answer the following questions: r How far is a von Kármán flow velocity time series from an ARMA(1, 0) model? r Can inhomogeneous anisotropic turbulence be better described by other ARMA(p, q) models and for which orders? r Is there a spatial organization of ARMA(p, q) reflecting the spatial distribution of velocity inhomogeneities? r Do different flow configurations correspond to different ARMA(p, q) models? r At a pure statistical level, is there an amount of information that the ARMA modeling can extract with respect to other techniques?

The main achievement of the paper is to suggest that these questions can be positively answered with a rather simple analysis. Moreover, once the order p, q is identified, one has immediately a criterion to build continuous stochastic models similar to the ones introduced in Refs. 11 and 22 for the quantities analyzed. Our aim is thus to define a general technique which can be then used to analyze and critically extract information from any turbulence experiment. In the present paper we underline the general procedure, leaving specific applications to future publications. The paper is organized as follows: first we give an overview to present the relevance of the Obukhov model, then describe ARMA(p, q) models for turbulence by giving a survey of their statistical and mathematical properties. Then we present the experimental set up and the quantities analyzed with the algorithm. Finally, we present and discuss the results obtained, outlining perspectives for the analysis of general turbulent datasets.

II. ARMA MODELS OF TURBULENCE

A. From the Obukhov model to an ARMA(1, 0) process

The celebrated phenomenological theories of Kolmogorov and Obukhov 1, 8 aimed to represent the complex phenomena of turbulence with a simple stochastic model. Thomson [START_REF] Thomson | Criteria for the selection of stochastic models of particle trajectories in turbulent flows[END_REF] was able to show that, in the inertial subrange, passive tracer Lagrangian velocities can be modeled by a Langevin equation (or Ornstein-Uhlenbeck process) with known coefficients; when discretizing this equation for simulation purposes, one can formally write it as an autoregressive process of order 1, usually denoted AR (1) or ARMA(1, 0). In particular, in the unidimensional case, the evolution of the velocity and of the position of a tracer particle (u, x) can be described by the stochastic differential equations:

du = a(x, u, t)dt + b(x, u, t)dW, (1)
dx = u dt, (2)
where dW are the increments of a Wiener process. In the same paper the determination of the coefficients a and b is discussed and it is found that, in Gaussian homogeneous turbulence, a = -u T L , where T L is the Lagrangian decorrelation timescale, while b = √ C 0 , where C 0 is a universal constant and is the mean kinetic energy dissipation rate. This can be written, as suggested in Ref. 24, in terms of macroscopic quantities:

= 2σ 2 u C 0 T L , ( 3
)
where σ 2 u is the fluid velocity variance (equal to the Eulerian variance) and can be seen as a measure of the turbulence intensity. Once the coefficients are known, one can write a discrete version of the Langevin equation ( 1): We are now considering a discrete-time stochastic difference equation, so we can use a discrete-time index t ∈ Z and, rearranging, Eq. ( 4) reads

u = - u T L t + C 0 W. ( 4
u t = 1 - t T L u t-1 + C 0 W. (5) Denoting 1 -t T L
= φ, √ C 0 t = σ and recalling that { W } are the increments of a Wiener process, the equation can be rewritten as follows:

u t = φu t-1 + ε t , ( 6
)
where {ε t } are independent variables, normally distributed. Equation ( 6) is the expression of an AR(1) process. To show that it is the simplest physical model which agrees both with Richardson law and the inertial range scaling proposed by Kolmogorov, it is sufficient to note that in an AR(1) process, the expected values of the velocity and the position scale in time, respectively, as

E[u 2 (t)] ∼ t, E[x 2 (t)] ∼ t 3 . ( 7
)
The second property is the Richardson law. Then, defining δu = E[u 2 (t)] and = E[x 2 (t)], we get from Eq. ( 7) δu ∼ 1/3 which can be seen as an equivalent of the Kolmogorov scaling.

B. Generalization: ARMA(p, q) model for turbulence

The ARMA(1, 0) leads to a Markovian evolution for the Lagrangian turbulent velocity, and is unable to describe the intermittency or memory that have been shown to exist in real flows. In most laboratory turbulent flows, available datasets are time series of values of a physical observable at a fixed point or obtained by tracking Lagrangian particles. In our case, time series are obtained at fixed points in space; in this work, no spatial velocity profiles are studied. This historically motivated the shift of paradigm from space velocity increments to time velocity increments defined as δu τ = u(t + τ )u(t) and motivated computations of the time structure function. Of course, in situations where measurements are made on the background of a strong mean velocity U, scale velocity increments and time velocity increments can be directly related through the Taylor hypothesis = Uτ . In situations such that the fluctuations are of the same order than the mean flow, however, the Taylor hypothesis fails. A suggestion has been made by Pinton and Labbé [START_REF] Pinton | Correction to the taylor hypothesis in swirling flows[END_REF] to then resort to a local Taylor Hypothesis, in which = dtu(t) where u is the local rms velocity. This is equivalent to consider a scale such that ∼ τ δu τ and may be seen as equivalent to modifying the space Kolmogorov refined hypothesis into a time hypothesis.

A natural generalization to take into account these features is thus to consider higher order ARMA(p, q) models, exhaustively treated, in example, in Ref. 26. A summary of useful notions about ARMA(p, q) modeling is provided in the Appendix. An ARMA(p, q) model corresponds to discrete time, stationary stochastic processes {X t } such that, for all t:

X t = p i=1 φ i X t-i + q j=1 ϑ j ε t-j + ε t . ( 8
)
Here {ε t } is assumed to be a white noise of variance σ 2 and the polynomials φ

(z) = 1 -φ 1 z t-1 -• • • -φ p z t-p and ϑ(z) = 1 + ϑ 1 z t-1 + • • • + ϑ q z t-q
, with z ∈ C, have no common factors. Notice the white noise assumption is a very general condition and X(t) will be normally distributed, resulting by a linear combination of independent and identically distributed random variables. From a physical point of view Eq. ( 8) is the natural extension of the ARMA(1,0) model corresponding to the Obukhov model by introducing a temporal memory structure: Intuitively, the autoregressive part of the process expresses a dependence of the value of the process at time t on a linear combination of its own p previous values, while the moving average component introduces, at time t, a linear dependence of the q previous values of the noise term. The quantification of memory effect in real turbulent flows will then be made through fits of the data by an ARMA(p, q) model, and measurements of how far this model is from the ARMA(1,0) model. For this, we first need to define the notion of best ARMA(p, q) fit, and then the notion of distance between the ARMA(1,0) and a given ARMA(p, q).

C. Model selection and characterization via correlation analysis

The main idea of time series analysis through ARMA models is to select the linear model that fits the data in the most parsimonious way, so that diagnosis of the nature of the generating process, forecasting or Monte Carlo simulations can be performed. Model selection is a non-trivial step of the procedure that can be addressed essentially by two means: through correlation analysis or through information based criteria, such as the Bayesian Information Criterion (BIC) described below.

For very simple processes, one can get access to the time dependence structure through computation of the auto-correlation function (ACF) and of the partial autocorrelation function (PACF) formally defined in the Appendix. In particular, for a MA(q) process, the theoretical ACF is characterized by q non-zero peaks, while the PACF decays exponentially or as a damped trigonometric function; for AR(p) processes the PACF is characterized by p non-zero peaks while the ACF decays exponentially. Hence, this fact allows to rule out or confirm the validity of an AR(p) or MA(q) hypothesis by a simple inspection of the ACF and PACF.

In the general ARMA(p, q) case, the simple correlation analysis described just above is not insightful. The model choice and the parameters estimation can be assessed by using the procedure introduced in Ref. 27, which also takes into account more complicated (such as integrated and seasonal) models:

1. preliminary analysis: the series is plotted in order to identify possible trends in mean and variance or periodic behaviors. Since here we deal with physically stationary processes, no trends are expected; 2. identification on the basis of the estimated ACF and PACF (or applying information criteria, such as the BIC); 3. estimation through maximum likelihood techniques; 4. diagnostic checking, that is, testing the estimated sequence for residual correlations (and normality or other distributive hypotheses, if required).

In the following analysis we perform the second step of the procedure fitting an ensemble of models with different (p, q) couples; we then choose the ARMA(p, q) model with the lowest total order p + q producing not correlated residuals. The serial independence of the residuals series is tested as described in the Appendix. As already mentioned, this phase could be based on the value of the BIC. In this case, the information criterion is computed for each model: the best fit is the minimum BIC after the steepest descent. The two methods provide the same results. First of all, we tested them on a synthetic time series of 10 5 values simulated from an ARMA(3,1), obtaining a correct estimation of the model with both methods. To ensure that this technique is stable also for shorter time series, in Fig. 1 we show a BIC profile as a function of p and q for one of the analyzed velocity samples, consisting of 600 observations: both the methods lead to the choice of an ARMA (7,7). All the PIV time series have length n = 600, while the typical LDV sample size is n ∼ 5 × 10 5 .

D. A measure of distance from Kolmogorov theory based on the Bayesian information criterion

It is useful to concentrate the information obtained by analysis in a single index. We want to obtain a measure of the distance of the selected ARMA(p, q) model from the ARMA(1, 0), namely the Thomson-Obukhov model.

For a given dataset, the relative quality of a statistical model can be measured by the BIC, defined as

B I C = -2 ln L + k[ln(n) + ln(2π )], ( 9
)
where L is the likelihood function for the investigated model. For an exhaustive definition of this quantity, see Ref. 28. Since the likelihood function is maximized when the correct model is found, while goes to zero in case of misspecification, its logarithm grows for well-specified models, while diverges to -∞ otherwise. Thus, the first term globally tends to become negative or to assume small values once the best model form is identified. On the other hand, the second term grows with the number of parameters times the sample size: so it serves as a penalization for the number of parameters, in order to avoid overfitting. In brief, when testing an ensemble of models for a certain dataset, the best one is identified by the minimum value of the BIC. For a Gaussian ARMA(p, q) model, it is expressed as follows:

B I C(n, σ 2 , p, q) = (n -p -q) ln n σ 2 n -p -q + n(1 + ln √ 2π ) + (10)
+ ( p + q) ln n t=1 X 2 t -n σ 2 p + q .
Notice that n is fixed by the experiment. The sample variance σ 2 is computed from the sample and is a series-specific quantity. Thus, in order to obtain a meaningful definition of the distance from Kolmogorov model, the BIC(n, σ 2 , p, q) must be normalized with respect to the Obukhov case BIC(n, σ 2 , 1, 0):

ϒ = 1 -exp |B I C(n, σ 2 , p + 1, q) -B I C(n, σ 2 , 1, 0)| /n, 0 ≤ ϒ ≤ 1. (11)
This quantity tends to zero if the dataset is well described by the Obukhov model and tends to one in the opposite case. We introduce the p + 1 correction to magnify small ϒ values.

III. EXPERIMENTAL SET-UP AND DATA PROCESSING

In order to illustrate and apply these concepts, we have worked with a specific axisymmetric turbulent flow: the von Kármán flow generated by two counter-rotating impellers in a cylindrical vessel. The experimental set-up is described in Refs. 21 and 29. Here, we consider a configuration where the disks are exactly counter-rotating at frequency f 1 = f 2 = F, in the two forcing conditions associated with the concave (resp., convex) face of the blades going forward, denoted in the sequel by sense (-) (resp., (+)). The resulting mean velocity fields are quite similar, with two toric recirculations separated by a mean shear layer (see Fig. 2). The forcing conditions however strongly influence the level of fluctuations, which are much higher in the (-) case. The working fluid is water, with viscosity ν = 1.0 × 10 -6 m 2 s -1 . The Reynolds number is defined as

Re = 2π F R 2 ν -1 ,
where R is the cylinder radius. We introduce a cylindrical system of coordinates x = (R, ϕ, Z ) with its origin at the center of the cylinder and the z-axis aligned with the impeller's rotation axis (see Fig. 1 In the sequel, we analyze both local and global observables. As local observable, we will consider the time series of the modulus of the velocity fields

V ( x, t) = [u( x, t), v( x, t), w( x, t)]
obtained by PIV measurements. Here u is the component in the PIV plane described in terms of R, the radial distance from the center of the cylinder; v in terms of Z, the vertical distance from the center, and w is the normal component to the PIV plane (the azimuthal velocity in cylindrical coordinates). For the comparisons between PIV and LDV measurements we will consider the normal component only w. We will also address two important aspects of statistical modeling of turbulence: the appearance of multifractal cascades, usually studied via the computation of the Hurst Exponents, and the role of phase randomization, which permit to isolate the effects of intermittency related only to the phase of the signals.

As global observables, we consider first the normalized kinetic energy introduced by Cortet et al.:

6 δ(t) = V 2 (t) V 2 . ( 12
)
Here the brackets indicate the spatial average, and the bar a time average. δ(t) represents the ratio of the total kinetic energy of the instantaneous flow to the total kinetic energy of the mean flow. As a second global observable, we also consider the torque time series C 1 (t) and C 2 (t) experienced by the two motors. The goal is to compare which part of information about the flow is carried by observables built using the velocity fields (such as δ(t)) and which is carried by dissipation measurements such as the torques.

The impeller speed F and the applied torques C 1 and C 2 are related to the average dissipation rate in the experiment, , through the injected power P. The typical kinetic energy in the experiment, V 2 , can be directly computed from the PIV data. Knowing these global quantities, it is possible to obtain rough estimates of two typical quantities of the turbulent flow, the Kolmogorov typical length and time scales η and t η and the Taylor typical scale λ, using the dimensional analysis inspired by Kolmogorov [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers[END_REF] and the identities of Taylor. [START_REF] Taylor | Statistical theory of turbulence[END_REF] These relations are only valid for homogeneous and isotropic turbulence, which is not the case here: we will use them anyway to get rough estimates, presented in [START_REF] Taylor | Statistical theory of turbulence[END_REF] Data are obtained with curved blades in both rotation senses, for f 1 = f 2 = 5 Hz, using water. Two turbulent states coexist in the (-) sense due to hysteresis, [START_REF] Ravelet | Multistability and memory effect in a highly turbulent flow: Experimental evidence for a global bifurcation[END_REF] one with one recirculation cell and another one with two cells. The one-cell state is unstable in the (+) direction for

f 1 = f 2 . Sense Cells (m 2 s -3 ) V 2 (m 2 s -2 ) η (m) t η (s) λ (m) ( -) 1

IV. RESULTS

We begin the analysis of the datasets by showing how the ARMA procedure, described in Sec. II, works on two velocity series extracted at two different locations from the same PIV experiment at Re = 10 5 and for the (+) sense of rotation. The series | V (t)| = V (t) [START_REF] Kato | Laser doppler velocimeter[END_REF] and their ACF and PACF are represented in Fig. 3. They have been obtained by sampling the data at 15 Hz. By analyzing the structure of the ACFs and the PACFs, one observes immediately that they do not consist of a small number of discrete peaks out of the confidence bands. This excludes the possibility that the series can be represented by pure AR(p) or MA(q) processes. Moreover, it is clear that a by-eye determination of the order (p, q) is not possible.

This result is consistent with the non-Markovian behavior used to describe the torque measurements via stochastic models in Refs. 11 and 22. By implementing our best fit procedure, we find that the best ARMA model to fit the data depends on the measurement points: the series on the left-hand side of Fig. 3 is fitted by an ARMA(1,1) model whereas the other one by an ARMA(4,2) model. This is of course not surprising, because the von Kármán flow is highly inhomogeneous. In the remaining of this section, we analyze the relationships between the flow inhomogeneous spatial structure and the ARMA fit structure by mapping the ARMA parameters.

A. Velocity fields

Let us now analyze the spatial structures obtained by applying the procedure described in Secs. II and III for a PIV field taken at Re = 2 × 10 5 , with mean velocity field provided in Fig. 2 for the (+) sense of rotation (left) and (-) (right). The two pictures look extremely similar: one can immediately recognize the cells structure of the flows described in Sec. III.

A full overview of the quantities computed by using the ARMA analysis is presented in Fig. 4 for the (+) sense of rotation and in Fig. 5 for the (-) sense. Obviously, even if the four cells structure presented in Fig. 2 (top left panels) and the standard deviation of | V (t)|(top right panels) show remarkable differences between the two configurations. In the (+) sense of rotation the four cells structure is appreciable whereas higher mean values and fluctuations of | V (t)| are recorded in the proximity of the wall of the cylinder in the (-) rotation. In this latter configuration, the fluid is pushed to the side of the cylinder and higher turbulent fluctuations are registered, as described in Ref. 5. Let us now analyze what happens to the quantities introduced by the ARMA analysis.

Order of ARMA

We start by describing the behavior of the total order O = p + q of the processes fitted for each time series (middle left panels). In both senses of rotation, the highest orders are concentrated near the impellers. However, differences appear in the other regions of the domain. In the (+) sense of rotation, the highest orders are found in correspondence to the highest fluctuations. Near the center of the cylinder the orders are low and, for some of the series, the signal is indistinguishable from noise (order 0). In the (-) set up, the highest level of turbulent fluctuations contribute to homogenize the behavior in larger areas such that a weak four cells structure is recognizable. This effect is probably linked to the presence of more homogeneous fluctuations in the flow (top right panel of Fig. 5). Even if the O color scale has been limited at p + q = 6 for comparison with the (+) situation, we that much higher orders appear in the (-) setups the impellers and, locally, at the walls of the cylinder. High p, q orders are directly connected to the vortices introduced by the rotations of the impellers and whose appearance is explainable in terms of Goertler instabilities. [START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF] We will see in Sec. IV D that these effects are recovered also for global observables.

Distance from Kolmogorov model

The difference between the (-) and (+) configuration is also highlighted by the results of the ϒ computations reported in the middle right panels of Figs. 4 and5. The distance from Kolmogorov model is lowest and almost zero near the boundaries, where the fluctuations are modest, and increases towards the center.

Evident differences appear if one compares ϒ values for (+) and (-) senses of rotation. As expected, the highest values are found in the (+) case, which is the one preserving a spatial four cells pattern in the fluctuations. This suggests that the coherent structures visualized using bubble air seedings are responsible for deviations from the Obukhov model.

In the (-) set-up, the region of values of ϒ ≥ 0.1 clearly traces the area with maximal azimuthal velocities. We have therefore a clear connection between coherent structures. In the present data set, we do not observe obvious signature of the influence of the shear layer dynamics. However, by using a much larger data set, we have been able to evidence the signature of the wandering of the shear layer in between to metastable position. This is reported in Ref. 33.

Physical interpretation of ARMA(p, q) coefficients

The bottom panels refer to the sum of the coefficients = p i=1 |φ i | (bottom left panels) and = q i=1 |ϑ i | (bottom right panels). and may be regarded as a representation of the total persistence of the phenomena, i.e., how much the system remembers of its past history. In order to get a better understanding of this idea and thus obtain a physical interpretation of the AR and MA parts of the process, we exploit, once again, Thomson's model. Equation (1) implies a Markovian evolution of the Lagrangian velocity in the inertial subrange, which is linked to an exponential behavior of the ACF:

ρ(t) ∼ e -t T L .
We can observe that the first term of the rhs of the equation contains some information about the global correlation structure of the process, which is even more evident in the discrete time (see Eq. ( 5)), since the autoregressive coefficient φ = (1 -t T L ) is the Taylor expansion of the exponential ACF. We have already mentioned that |φ| is a measure of the persistence of the process; here, this persistence is driven by the large eddies, since T L is the Lagrangian decorrelation time scale, which is T L t, with t lying in the inertial subrange. On the other hand, since the noise term in Eq. ( 1) is a standard Brownian motion, the innovations of the time series in Eq. ( 5) W , are normally identically distributed and δ-correlated. Thus, the second term of the equation is simply a noise driving the process, with no linear dependence between any couple of values ( W i , W j ) with i = j; this assumption is due to the fact that the stochastic kicks come from the viscous eddies, which live in the viscous time scale τ η t: since the particle samples turbulent field with a frequency linked to a characteristic time ∼ t, the viscous eddies are completely uncorrelated between two steps. This means that here q = 0 and the innovations of the process are a pure (Gaussian) noise, and all the information is contained in the mean and the variance.

In general, we may say that the AR(p) part of an ARMA(p, q) process is linked to the contribution of the large scales and represents the persistence of the process. Notice that, if p = 1, it must be |φ| < 1 in order to satisfy the stationarity condition; if the process is more persistent than an AR(1) with |φ| < 1, higher values of p are required to explain all the correlation coming from the large scales. Analogous considerations hold for the MA(q) part: if q > 0, a linear combination of previous values of the noise appears in the equation, introducing a correlation structure in the innovation term, i.e., a higher persistence of the noisy contributions. This means that the small eddies do not decorrelate completely between two sampling times, so we should assume to have eddies at all scales.

Comparison with high-resolution datasets (LDV)

The results obtained with the PIV technique must be validated and checked against higher temporal resolution datasets. In fact, although the possibility of defining a distance from the Kolmogorov model in the physical space rather than in the Fourier space seems appealing, we must be sure that the results obtained with the ARMA analysis are stable with respect to an increase of resolution.

In the previous discussion, we have pointed out that the Obukhov model, representing homogeneous and isotropic turbulence, can be written as an ARMA(1,0) model. This corresponds to have a purely power-law spectrum which does not contain any other features than the decay predicted by Kolmogorov. Since Eq. ( 5) contains an explicit dependence on t only for the coefficient φ, we do not expect to see a change in the order of the process when increasing the resolution, but rather changes in φ and ϑ coefficients. This consideration holds unless the spectrum changes slope or has peaks for some of the frequencies we add to the spectrum by increasing the resolution. In this case we expect to see also a change of the autoregressive and moving average polynomials.

In our analysis, we compared the PIV data for the (-) sense of rotation with the LDV experiment performed in the same conditions. Since for the LDV series only the w component is measured, we will compare this quantity to the same recorded for the PIV experiments. The LDV data allow for exploring frequencies of order of the kHz, whereas the PIV is limited to a frequency of 15 Hz, so that we extend significantly the range of frequencies analyzed. Whereas the time resolution of the LDV data is very high, the spatial resolution is indeed low: we have w measurements only at the 18 points represented by the red crosses in the top-left panel of Fig. 6. For this reason, the quantities obtained from the LDV analysis (left panels of Fig. 6) have been interpolated on a finer spatial grid. Anyway, the level of details remains lower if compared with the PIV results (right panels of Fig. 6). The top panels of Fig. 6 show a comparison between the averaged velocity field as obtained from the LDV and the PIV analysis. They both show not only the familiar cells structure, but also that the order of magnitude of the velocity fields is extremely close for the two different techniques. The analysis of the orders O (reported in the central panels of Fig. 6) shows consistency between the two techniques: the highest orders are located at the walls of the cylinder and near the impellers. Moreover, if we average the total order on all the available points, we get O = 2.4 ± 0.7 for the LDV and O = 1.8 ± 0.8, values compatible within a standard deviation. Finally, the analysis of ϒ (lower panels of Fig. 6) reveals that the maxima are located, both for the PIV and the LDV, near the walls of the cylinder around Z = 0. By computing the average of ϒ over all the points we find 0.03 ± 0.02 for the LDV data and 0.02 ± 0.02 for the PIV, again consistent within a standard deviation.

B. ARMA analysis of phase randomized data: Phase intermittency

Phase randomization is often used in turbulence to destroy the intermittency effects related to Fourier phases while preserving the intermittency effect related to Fourier amplitude: by applying such procedure one preserves up to the second order statistics (covariance and spectrum). It is therefore interesting to apply the ARMA analysis to phase randomized data sets to quantify the relative influence of phase and amplitude intermittency in turbulence.

A simple and efficient way to perform this phase randomization is to compute the Fourier transform of the time series, then randomize the phase (while preserving the anti-symmetry of the phase with respect to the frequency variable resulting from the real nature of the data) and going back to the physical space by means of an inverse Fourier transform. In order to perform this task, we have used the MATLAB code provided by Carlos Gias, based on the procedure described in Ref. 35.

The results we present correspond to the (-) sense of rotation for the PIV data already analyzed in Sec. IV. After generating surrogate velocity data, we compute ϒ for the phase-randomized data and compare it with the original one. This is done in the lower panel of Fig. 7, with the original ϒ (left) and the ϒ for phase-randomized data (right). Both panels show the same structure meaning that most of the contribution to the intermittency parameter is associated to intermittency amplitude through first and second order statistics (presumably through the advection and shearing effect of the large scale flow). To get information about phase intermittency, we subtract the intermittency index from the phase-randomized data to the original ones, obtaining ϒ. Results of such a difference are reported in the upper panel of Fig. intermittency. Its resulting spatial structure is highly organized, showing some association with vortices. This feature may be connected with the observation that anomalous scaling, in linearly advected hydrodynamical models, is connected to the existence of statistically preserved structures with highly complex geometrical properties. [START_REF] Biferale | Anomalous and dimensional scaling in anisotropic turbulence[END_REF][START_REF] Celani | Statistical geometry in scalar turbulence[END_REF] We leave this for future investigation.

C. Hurst exponents

A generalized version of the first equation in the system (1) can be written as

du = a(x, u, t)dt + b(x, u, t)dW 2H , ( 13
)
where dW H is the increment of a fractional Brownian motion (fBm) and H is the so-called Hurst exponent. [START_REF] Beran | Statistics for Long-Memory Processes[END_REF] The fBm, first introduced by Mandelbrot and Van Ness, [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF] is a generalization of Brownian motion where the increments are not independent. It has zero mean and the following covariance function:

E[W H (t)W H (s)] = 1 2 (|t| 2H + |s| 2H -|t -s| 2H ).
The exponent H is a real number in (0, 1) and its value determines the memory of the stochastic process. For H = 1/2, the standard Brownian motion is recovered. For 0 < H < 1/2 the process is anti-persistent, i.e., an increase will most likely be followed by a decrease or vice-versa. Finally, for 1/2 < H < 1, the series is persistent, i.e., increases generally follow increases. We want to investigate if the behavior displayed by the total order and total persistence of the ARMA(p, q) and by the distance index ϒ can be better explained by the fractional nature of the underlying stochastic process. In order to do this, we compute H for the | V (t)| data in the (+) sense of rotation and we compare it to ϒ in Fig. 8. The computation of the Hurst exponents follow the methods presented in Ref. 40 which we found to be all consistent with each other. Since ϒ values span 5 orders of magnitude while H is always of order 1, we consider the log 10 (ϒ). Not only the spatial structure of H and log 10 (ϒ) are very similar (upper panels of Fig. 8), but also a linear relation can be found between these two quantities (lower panel of the same figure). The linear correlation coefficient is r = 0.70 and these results hold also for the LDV experiments with almost identical fit coefficients and r = 0.81. From this analysis we argue that a fBm description of the phenomenon might be used to explain the nature of the correlations in the series and it could be useful to further improve the modeling of inhomogeneous and anisotropic turbulence. However, the results show that ϒ, based only on an ARMA(p, q) estimation, is equally effective in quantifying deviations from the Kolmogorov model which could not be due to fBm effects.

D. Global observables

An interesting question to address when dealing with spatial-temporal extended systems, is how the information on the single trajectories is transmitted to integrated quantities. In particular, one may ask whether the differences found in the ARMA analysis for the local observables of PIV fields are preserved for scalar quantities, i.e., if high ARMA(p, q) orders found locally in the proximity of the impellers and the cylinder walls give a contribution to global observables or whether they average out. In this section we present results obtained for the quantity δ(t) introduced in Eq. (12). We have further tested that our results are independent of the choice of the global observable, whether derived from PIV measurements -Angular momentum -or measured independently like torques measurements.

We have carried out the analysis on global observables at several Reynolds numbers around Re = 10 5 , that is in a fully turbulent regime. The typical behavior of the time series of δ(t) is represented in Fig. 9 for the (+) sense of rotation (left), and the (-) one (right). The top panel refers to the time series obtained by averaging the spatial velocity fields and shows no particular differences at first sight, as we have seen for the examples of the velocity series shown in Fig. 3. However, the ACF (middle panels) and PACF (lower panels) are remarkably different. The ACF of the (+) sense of rotation decays quickly and the PACF shows only one peak significantly different from zero: an ARMA(1,0) is enough to explain the correlation structure. On the other hand, an oscillatory behavior of both the ACF and PACF is clearly recognizable for the (-) rotation. The orders p, q needed to decorrelate the latter time series are higher in the (-) rotation, namely p = 2, q = 1. These results hold generally by varying Re and changing observables and point to the intrinsic differences between the two senses of rotation.

One can notice that some characteristic features appear in the ACF of the (-) sense of rotation. We can speculate that, for scalar quantities, the highest orders get averaged out if their contribution is substantially different at different (r, z) as it happens for the (+) rotation. However, when the same kind of features are present in the ACF and PACF for series at different (r, z) the contribution sums up and is well visible in the behavior of global observables.

V. MULTISTABILITY

Another interesting question is whether the application of ARMA techniques to turbulence is helpful to discriminate between different stability regimes. The simple guess is that by increasing the instability of a configuration, higher orders O arise as we introduce in the system new time scales linked to the presence of nearby attracting states. In order to check this idea, let us consider again a von Kármán swirling flow with the same geometry described before, with the Reynolds number fixed at Re ∼ 10 5 . In this system, one can impose either the speeds f 1 and f 2 of the motors or the torques C 1 , C 2 and define two natural dimensionless quantities:

θ = ( f 1 -f 2 )/( f 1 + f 2 ), γ = (C 1 -C 2 )/(C 1 + C 2 )
which are, respectively, the reduced impeller speed difference and the reduced shaft torque difference. In Ref. 21, the authors found that different forcing conditions change the nature of the stability of the steady states. Here we complete the results represented in Fig. 1

Speed control:

In this case all the turbulent flows are steady. By plotting the averaged θ 's and γ 's measured in several experiments, one obtains the diagram shown in Fig. 10. The colors refer to different O (Fig. 10, upper panel), and R (Fig. 10, lower panel). Starting both impellers at θ 0 leads to a marginally stable state, which consists of two symmetric recirculation cells separated by a shear layer. If one waits enough time, a fluctuation may force a jump of the system to one of the two bifurcated states represented by the red points. The instability of the symmetric state is reflected by the order of the ARMA processes fitted for the series of γ (t) at γ θ 0. For this experiment we found O = 4 and R 3, values definitely larger than the ones found in the bifurcated states where always O = 1 R < 0.5. From the available data one can argue that the potential barrier-the repellor in dynamical system-at θ = ±0.1 is somehow impenetrable as we do not get any increase in O and R for the series γ (t) recorded at such values of θ .

Torque control: By imposing the torque control one gets access to new attracting states, located in correspondence to the repellor found in the speed control. The results for the torque control have been obtained by analyzing time series of θ (t) and the results in terms of O and R are reported in Fig. 11. Before commenting on the new states, we begin by analyzing the states which are attracting in both the configurations. In Ref. 21, the authors assert that the properties of the attracting states in the speed control set up are analogous to the ones found for the torque control. However, by applying the ARMA analysis, we found, as one would expect, remarkable differences. In particular, the symmetric state (which was marginally stable in the speed control) is now stable as one can go to the bifurcated states in a continuous way. This is confirmed by the low values of O and R found for the torque control where O = 1 (Fig. 11(a), found in correspondence of γ 0.005 and θ 0) and R < 1 (Fig. 11(b)). The bifurcated states have a different characteristic order (typically O = 2) and persistence (typically 1 < R < 2). These are not linked to the presence of transitions, as they persist further away from the unstable range of parameters. They are linked to the modifications in the dynamics induced by the change of control which affects also the stable regions in a fine way, which has not been discussed in Ref. 21 but is evident by applying the ARMA technique capabilities. Let us now comment on the new states which appear in the torque control in correspondence to the repellor for the speed control. These states feature multistability as detailed in Ref. 21. In terms of ARMA analysis they are characterized by higher orders (green branches in the upper panel of Fig. 11 with O = 3), and persistence (black branches of the lower panel in Fig. 11). Even if an order is found for multistable time series, an ARMA model cannot just be stuck to the data in this case, as it will not reproduce a multistable behavior, but rather a process with one stable state whose correlation properties are similar to the ones found for the multistable time series.

θ γ R ≤ 0.5 0.5 ≤ R ≤ 1 1 ≤ R ≤ 2 2 ≤ R ≤ 3 R ≥ 3 (a) (b)
This example clearly shows that one can find, far from the bifurcation, a typical ARMA(p, q) process (O = 1 for the speed control and O = 2 for the torque control) which describe the data-sets whereas p and q are evidently modified by the stability properties. We remark that one must pay extraordinary care when the goal is to find a model for the data-set.

VI. FINAL REMARKS

In this paper we have shown how to extract information from turbulent data-sets by applying an ARMA statistical analysis. Such analysis goes well beyond the analysis of the mean flow and of the fluctuations; in fact, it is possible to link the behavior of the recorded time series to a discrete version of a stochastic differential equation which is able to describe the correlation structure in the data-set. We have tested the method on data-sets produced by the experiments on the von Kármán swirling flow, already analyzed in several publications.

We have shown how the anisotropies and inhomogeneities present in real experiments as well as the finite resolution of the data-sets influence the order p, q of the process which better describes the data. We find that data are suitably described by ARMA(p, q) processes whose orders are generally different from the Obukhov model although with a very limited number of auto-regressive and moving average terms (generally p, q = 1 or 2). We have introduced a new index 0 ≤ ϒ ≤ 1 to measure and quantify this difference. The value of ϒ increases in areas where large scale coherent structures are present. It would be interesting to rely the statistics of ϒ to the computation of refined statistics of velocity increments, possible only for time high-resolved experiments. [START_REF] Arneodo | Structure functions in turbulence, in various flow configurations, at Reynolds number between 30 and 5000, using extended self-similarity[END_REF] In particular, we aim to compare ϒ with classical intermittency parameters based on structure functions. This idea follows from the hypothesis first proposed by Laval et al. [START_REF] Laval | Nonlocality and intermittency in three-dimensional turbulence[END_REF] that intermittency propagates in direct interactions between large and small scales, rather than in cascades. Preliminary analysis carried on the LDV data-sets show that there is a linear proportionality between ϒ and the classical intermittency parameters defined on the velocity increments. The inhomogeneous structure of the PIV experiments is reflected by the range of different orders p, q found in our analysis: a great part of the flow can be described in terms of noise, whereas higher orders concentrate around cells ((+) sense of rotation) or near the walls ((-) rotation). The analysis of global observables shows that most of the information about the local structure of the flow is preserved, including the differences found between the two senses of rotation. This correspondence between local and global quantities is very important and it will be further exploited for challenging systems for which only global observables are available as in the SHREK data-sets experiment with super-fluid Helium. [START_REF] Saint-Michel | Probing quantum and classical turbulence analogy through global bifurcations in a von Kármán liquid Helium experiment[END_REF] We have also checked our results against the change in time resolution by comparing them against the LDV experiments, whose average sampling frequency is order of the kHz. Not only the values of O found with this technique are consistent, but also the values of ϒ and the spatial structures observed. Finally, we have commented on the effects of multistability on the ARMA analysis by considering two different kind of forcings for the von Kármán experiment.

In the Obukhov model, the coefficients of the stochastic model are given from the turbulence theory, resulting in a simple Langevin equation which describes the process in the continuous time. Here, we have applied estimation methods to obtain a parametric description in the discrete time, but the passage to continuous time stochastic differential equations is not trivial for a general ARMA(p, q) process. Obtaining an expression of the model's parameters in terms of physical quantities of turbulence theory is presently not possible. In fact, it is likely that the MA(q) part of our processes represents the contribution of the shortest time scales detectable with the available techniques. One way to test this idea is to verify that only the structure of the MA(q) part of the process changes by changing only small scales feature of the flow. In order to do that, we are currently testing impellers with a fractal structure. Preliminary analysis show that MA(q) orders are different between fractal and non-fractal impellers whereas the AR(p) do not change. Details will be reported in a future publication.

Several generalizations of ARMA models exist and they allow for taking into account the possible multi-fractal behavior of turbulence. The comparison of the Hurst exponent and the ϒ index suggests that it will be interesting to extend the analysis to fractional integrated ARMA or ARFIMA(p, H, q) models. The first ones can be appropriated for studying problems of non-stationary turbulence [START_REF] Connaughton | Non-stationary spectra of local wave turbulence[END_REF] whereas a SARMA models analysis could be suitable for studying problems of wave turbulence. [START_REF] Nazarenko | Wave Turbulence[END_REF]

(Received 25

 25 March 2014; accepted 16 September 2014; published online 6 October 2014)

FIG. 1 . 6 Faranda

 16 FIG.1. BIC values resulting from different fits of ARMA(p,q) model for a velocity time series consisting of 600 observations.

105101- 7 Faranda

 7 FIG. 2. Structure of the mean velocity field for θ = 0. The arrows represents directions and intensities of the velocity components in the PIV plan averaged over time ( ū, v). The orthogonal component w is represented by the color scale. Left: (+) sense of rotation. Right: (-) sense of rotation.

FIG. 3 .

 3 FIG.3. Two time series of |V (t)| (upper panels) with their respective ACF functions (middle panels) and PACF (lower panels). Re 10 5 , (+) sense of rotation. Blue lines in the ACFs and PACFs represent the confidence bands at the 95% confidence level. Sample frequency: 15 Hz. X-axis is in sample index.

  FIG. 5. ARMA analysis for the (-) sense of rotation. Top left: | V |. Top right: | V (t)| standard deviation. Center left: Total order O found by fitting an ARMA(p,q) to the | V (t)| data. Center right: ϒ for the | V (t)| data. Bottom left: Sum of the autoregressive coefficient . Bottom right: Sum of the moving average coefficients Theta.

  FIG. 6. Comparison between the ARMA analysis for LDV data (left) and PIV data (right) for the azimuthal component of the velocity w. w (upper panels), total order O (central panels) and distance from Kolmogorov model ϒ (lower panels); (-) sense of rotation; Re 10 5 . The red crosses in the top-left panel show the locations of the measurement points for the LDV experiment.

1 FIG. 7 .

 17 FIG. 7. Upper panel: Difference between ϒ computed for the | V (t)| data for the (-) sense of rotation and ϒ on surrogate date of the same experiment obtained with a phase randomization procedure. Lower panel: ϒ for the original data (left) and for the surrogates (right).

FIG. 8 .

 8 FIG.8. Upper panels: Comparison between the Hurst exponents H (left) and the deviation from the Kolmogorov model in log-scale log 10 (ϒ) (right) found by fitting an ARMA(p,q) to the | V (t)| data for the (+) sense of rotation. Lower panel: Scatter plot of the Hurst Exponent H and log 10 (ϒ). The red solid line shows a linear fit to the data.

FIG. 9 .

 9 FIG. 9. Time series of δ(t) (upper panels) with their respective ACF functions (middle panels) and PACF (lower panels). Left: (+) sense of rotation. Right: (-) sense of rotation. Re 10 5 . Blue lines in the ACFs and PACFs represent the confidence bands at the 95% confidence level. Sample frequency: 15 Hz. X-axis is in sample index.

2 FIG. 10 .

 210 FIG.10. Total order O and total persistence R for the von Kármán experiment under the speed control. The points represent the averaged γ and θ obtained for each experiment.

FIG. 11 .

 11 FIG. 11. (a) Total order O and (b) total persistence R for the von Kármán experiment under the torque control. The points represent the averaged γ and θ obtained for each experiment.

  )

	105101-4	Faranda et al.	Phys. Fluids 26, 105101 (2014)

  of Ref. 30).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.49.25.36 On: Mon, 06 Oct 2014 12:04:43

TABLE I .

 I Typical length scales of the flow, using only global average quantities and usual turbulent identities.

  This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.49.25.36 On: Mon, 06 Oct 2014 12:04:43

	105101-17	Faranda et al.							Phys. Fluids 26, 105101 (2014)
		0.15
			O=1
		0.1	O=2
		0.05	O=3
		γ	O=4
		0
		-0.05
		-0.4 -0.1	-0.3	-0.2	-0.1	θ	0	0.1	0.2	0.3
		0.15
		0.1
		0.05
		0
		-0.05
		-0.4 -0.1	-0.3	-0.2	-0.1		0	0.1	0.2	0.3

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 193.49.25.36 On: Mon, 06 Oct 2014 12:04:43

ACKNOWLEDGMENTS

We acknowledge the other members of the VKE collaboration who performed the experiments: C. Wiertel and V. Padilla. We acknowledge two anonymous referees whose comments greatly improved the results and the consistency of this paper. D.F. acknowledges the support of a CNRS postdoctoral grant.

APPENDIX: THEORETICAL SURVEY ON ARMA MODELING

We have already mentioned, in the definition of the ARMA(p, q) process, that {X t } must be stationary. Usually, two definitions of stationarity are given when treating discrete time stochastic processes: strong stationarity implies stationarity of the whole joint probability distribution of the stochastic process, while weak stationarity requires the first two moments of the process to be finite and constant in time. The results about ARMA(p, q) processes are usually proved requiring weak stationarity, that is of course implied by strong stationarity; in our data analysis we will assume stationarity on a physical basis, by studying the system when the dynamics has reached well identifiable stationary states.

First of all, we observe that Eq. ( 8) can be written in a very compact form, introducing the backward operator B such that B j X t = X t-j , j ∈ Z:

If ϑ(z) ≡ 1 the process reduces to an AR(p), if φ(z) ≡ 1 the process is said to be a MA(q) process.

Two other important features of a discrete time stochastic process are causality and invertibility. Causality refers to the possibility of recovering the value of the process at time t as a function of the innovations ε s , with s ≤ t. Formally, {X t } is causal if there exists a succession of absolutely summable coefficients {ψ j } so that the process at time t can be written as

which implies, in terms of the auto-regressive polynomial:

Invertibility could be regarded as the property specular to causality, so {X t } is invertible if the series of the innovations {ε t } can be recovered from the process. This requires the existence of a succession of summable coefficients {π j }, which allows us to write

and in this case, condition (A4) is required on the moving-average polynomial:

In case of the presence of d unit roots in the auto-regressive polynomial, the process becomes non-stationary; however, the d-differenced process (1 -B) d X t can be stationary. In particular, if (1 -B) d X t is an ARMA(p, q) process, X t is said to be an ARIMA(p, d, q) process (where the "I" stands for integrated); the particular case of an AR(1) with φ = 1 reduces to the well-known random walk. Taking the differences of a time series is a drastic operation and a careful testing for the presence of unit roots must be performed if this kind of non-stationarity (also called stochastic trend) is supposed to exist. The most used test to this purpose is the Augmented Dickey-Fuller test; notice that a unit root in the moving average can also be taken into account through the hypothesis that the innovations are an integrated autoregressive process. Extensions to more complicated models can be found in literature, but these basics ARIMA processes are sufficient for the data analysis proposed in the present work.

The main idea of time series analysis through ARMA models is to select the linear model that fits the data in the most parsimonious way, so that diagnosis of the nature of the generating process, forecasting or Monte Carlo simulations can be performed. Model selection is a non-trivial step of the procedure and should be discussed after the introduction of some fundamental tools for the investigation of the time-dependence structure of the stochastic process. This issue can be addressed through spectral analysis, decomposition of the time series in trend, cyclical, periodical, and irregular components and, most of all, correlation analysis. In this approach the dependence structure is studied analyzing the global and the partial autocorrelation functions.

The (global) auto-covariance function (ACVF) at lag h of a zero-mean stochastic process is defined as

and, when normalized over the variance, gives the (global) ACF: The concept of PACF is less intuitive; formally, it can be written as

where φ hh is the last component of

and γ h = [γ (1), . . . , γ (h)] . In practice, the PACF quantifies the correlation between the prediction errors at lag 0 and h, given that it can be shown that the conditional expected value is the best linear predictor:

For very simple processes, the ACF and the PACF give strong hints about the time dependence structure. In particular, for a MA(q) process, the theoretical ACF is characterized by q non-zero peaks, while the PACF decays exponentially or as a damped trigonometric function; for AR(p) processes specular considerations are valid. On the (+), if the process is characterized by both autoregressive and moving-average polynomials, few information can be obtained by a simple byeye evaluation of the ACF and PACF. In this case, statistical information criteria must be used; one of the most widely known is the BIC, which will be presented after the introduction of the concept of estimation for ARMA(p, q) models.

Given a parametric hypothesis ARMA(p, q) for a time series, the corresponding discretetime equation is fitted to the data and all the parameters are estimated with maximum likelihood techniques, well described also in the more practical volume. [START_REF] Brockwell | Introduction to Time Series and Forecasting[END_REF] At this point, two sets { φ j } p j=1 and { θ j } q j=1 of estimated parameters are available, as well as a time series of estimated residuals {ε} of the same length of the original time series; if the tested ARMA(p, q) fits the data, {ε} must be a sequence of independent random variables. Notice that, if the orders p and q are too high, the time series is over-fitted, so the analyst must be careful in choosing the most parsimonious model in terms of number of parameters. Thus, once the residual sequence has been obtained, inference must be made on the null hypothesis H 0 of uncorrelated residuals. At this point, the sample ACF ρ( j) is computed; then, one of the most used test statistics is the Ljung-Box Test:

where n is the length of the time series and h is a fixed number of lags at which the sample ACF is computed. If H 0 is not rejected at a given level (usually α = 0.01 or α = 0.05), the tested ARMA(p, q) fits the time series. As already mentioned, in case of complex or high-order processes, ACF and PACF are not sufficient to obtain a hint on the possible order (p, q) at glance; in this case, some different hypothetical ARMA(p, q) models can be fitted and for each one the BIC is computed:

This quantity is minimized by the most parsimonious model providing a good fit to the time series. However, an other possibility is to fit many different ARMA(p, q) models and choose the one for which the null hypothesis of uncorrelated residuals is not rejected and the total number (p + q) of parameters is minimum.


Christophe Baudet

Michel Bon-Mardion

Patrick Bonnay

A Braslau

Bernard Castaing

Francesca Chillà

Laurent Chevillard

François Daviaud

Pantxo Diribarne

Bérengère Dubrulle

Lo- Cal Velocity

M Bon Mardion

D Durì

D Faranda

B Gallet

M Gibert

A Girard

B Hébral

T Lehner

I Moukharski

J M Poncet

J.-P Moro

P.-E Roche

B Rousset

E Rusaouën

B Saint-Michel

J Salort

E-W Saw

K Steiros

C Wiertel-Gasquet

Local velocity measurements in the Shrek experiment at high reynolds number

The SHREK experiment (SuperFluid von Kármán experiment) has been designed to compare properties of superfluid turbulence with ordinary turbulence. The SHREK experiment (see Fig. 1) makes particularly use of the phase transition experienced by helium at T λ ≈ 2.2 K, thus allowing to make experiments at very high Reynolds values with normal (above T λ ) and superfluid helium (below T λ ) [START_REF] Maurer | Local investigation of superfluid turbulence[END_REF][START_REF] Salort | Turbulent velocity spectra in superfluid flows[END_REF]. The experimental setup is described in [START_REF] Rousset | Superfluid high reynolds von kármán experiment[END_REF]. It consists of a cylinder of radius R = 0.39 m and height H = 1.2 m (Fig. 1). The fluid is mechanically stirred by a pair of coaxial impellers. The impellers are disks of radius 0.925R, fitted with 8 radial blades of height 0.2R and curvature radius 0.4625R. The impellers rotation rate f 1 and f 2 can be varied independently from 0.1 to 2 Hz and the total dissipated power can range from 120W at 1.6 K to 400W above 1.9 K. This experiment benefits from high flexibility of flow conditions, due to the large variation of helium properties over the available temperature range (1.6 K to 5 K). Both superfluid and normal turbulence measurements are possible in the same experiment, with adjustable fraction of the superfluid component (from ≈ 85% at 1.6 K to 0% above T λ ). Acquisitions performed in HeJet experiment [START_REF] Duri | Liquid helium inertial jet for comparative study of classical and quantum turbulence[END_REF] First we bring new insight in the understanding of the signal of hot-wires in He II (see [START_REF] Duri | Heat transfer enhancement on thin wires in superfluid helium forced flows[END_REF]). In particular, we show that the inter-correlation coefficient of the signal of a hot-wire with the signal of a Pitot tube (for which the design can be found in [START_REF] Salort | Turbulent velocity spectra in superfluid flows[END_REF]), is the same in He II as in He I where the signal can reliably be interpreted as a velocity (see Fig. 2). We also show preliminary velocity measurements obtained in SHREK using this kind of hot-wire. Second, we report the first results obtained using novel miniature total head pressure tube (Pitot tube), obtained for a large range of rotation frequencies f 1 and f 2 , in both normal and superfluid regime. The Pitot Helmholtz frequency is ≈450 Hz giving useful bandwidth somewhat higher. The sensitivity is in the 10 -7 m 2 /s 2 range and can be further improved by combination of 2-amplifiers. The results are compared with results obtained in a scale 1:4 experiment, operated in Saclay with water [START_REF] Saint-Michel | Influence of reynolds number and forcing type in a turbulent von kàrmàn flow[END_REF].

Figure 1 .Figure 2 .

 12 Figure 1. Schematic view of the experimental setup and the impellers blade profile

Acknowledgment

This work was supported by the French National Research Agency grant ANR-09-BLAN-0094-01 and by the European Community Framework Programme 7, EuHIT -European High-performance Infrastructures in Turbulence, grant agreement no. 312778.


Marcelo M Wanderley
email: marcelo.wanderley@mcgill.ca

Joseph Malloch
email: joseph.malloch@inria.fr

Jérémie Garcia
email: j.garcia@gold.ac.uk

Wendy E Mackay
email: mackay@lri.fr

Michel Beaudouin-Lafon

Stéphane Huot
email: stephane.huot@inria.fr

Human Computer Interaction meets Computer Music: The MIDWAY Project

Keywords: Musical Interaction, Human-Computer Interaction, Digital Musical Instruments, Instrumental Interaction. ACM Classification H.5.2. User interfaces (Interaction styles); H.5.5. Sound and music computing (methodologies and techniques);

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Although Human-Computer Interaction and Music both rely heavily on interactive technologies, their goals sometimes differ. For example, HCI often focuses on improving human performance with interactive technologies, measured in terms of efficiency and accuracy. Music creation lacks such easily quantifiable and measurable goals, and instead focuses on concepts such as creativity, engagement, personalization, and appropriation, all difficult to account for in design.

We believe that the creative context of music provides opportunities for putting cutting-edge HCI models and tools into practice. Expert musicians push the boundaries of system design through their personalization and appropriation of music applications, rendering standard HCI performance measures insufficient for evaluating musical tools. We argue that computer music technology should support and embed existing design models and methodologies, and that new designs built from a well-defined design space will, in turn, facilitate validation and evaluation, as well as reuse and appropriation from other fields, making exploration and extrapolation possible [START_REF] Huot | Designeering Interaction': A Missing Link in the Evolution of Human-Computer Interaction[END_REF].

Methods and Models

HCI researchers have proposed many models to support the design and evaluation of interactive systems. Two are particularly relevant for musical interaction, with direct implications for our work:

Rasmussen's Human Information Processing [START_REF] Malloch | Towards a new conceptual framework for digital musical instruments[END_REF][START_REF] Rasmussen | Information processing and humanmachine interaction: an approach to cognitive engineering[END_REF] is a useful theoretical framework for making sense of the various interaction possibilities in music (Fig. 1). It has implications for the choice of mappings, e.g., from the static and deterministic mappings of digital musical instruments (DMIs) to the dynamic mappings of interactive music systems. It also suggests how to design feedback for performers, including continuous signals for skill-based performance, signs for rule-based interaction, and symbols for model-based behavior.

Beaudouin-Lafon's Instrumental Interaction [START_REF] Beaudouin-Lafon | Instrumental inter-action: an interaction model for designing post-WIMP user interfaces[END_REF] describes how users interact with objects of interest, mediated by 'Interaction Instruments' similar to interaction with physical tools. Instruments are treated as first class objects, and can be defined in terms of reification, polymorphism and reuse [START_REF] Beaudouin-Lafon | Polymorphism and Reuse: Three Principles for Designing Visual Interfaces[END_REF]. Instrumental interaction relates more directly to musical interaction than earlier HCI models such as direct manipulation.

Even so, musical instruments are usually inflexible: DMIs typically involve static mappings between sensor outputs and sound synthesis inputs that performers must adapt to. Although this limits their use, it allows musicians to master their instruments, sometimes with virtuosic skill. However, musicians also co-adapt [START_REF] Mackay | Responding to cognitive overload: Coadaptation between users and technology[END_REF]not only do they adapt their behavior to the instrument as they learn it, they also adapt or modify it to achieve creative goals. We hope to create co-adaptive instruments that are explicitly designed for appropriation by musicians, using rule-based approaches together with the principles of polymorphism and reuse.

Examples of Applications

We describe two projects for which we designed modular and reusable tools that meet the idiosyncratic needs of musicians. These applications illustrate and support concepts defined in the above models, opening new possibilities for flexible musical applications.

LibMapper & ICon

Input Configurator (ICon) [START_REF] Appert | FlowStates: prototypage d'applications interactives avec des flots de données et des machines à états[END_REF] and libmapper [START_REF] Malloch | Distributed tools for interactive design of heterogeneous signal networks[END_REF] are open-source software tools intended for the design and (re)configuration of modular interactive systems. Both focus on mapping and configuration as a top-level task ICon is an interactive data-flow editor in which various input and interaction resources, such as input devices and communication protocols, can be instantiated as data-flow processing devices, and connected to each other and to interactive applications. Libmapper was designed to support the creation of DMIs. Applications and input devices declare their local resources, which are remotely discovered and connected over a local network. Various session-management tools can be used to interact with the resulting distributed network to add, modify or remove connections between producers and consumers of real-time data.

ICon and libmapper are complementary. ICon features a much richer visual programming interface with a large, extensible library of data processing devices, as well as both data-flow and state-machine approaches for describing and prototyping advanced interaction techniques [START_REF] Appert | FlowStates: prototypage d'applications interactives avec des flots de données et des machines à états[END_REF]. It follows the principles of Instrumental Interaction by reifying interaction techniques into dataflow processing devices and configurations that can be manipulated as first class objects and applied to other contexts or systems (polymorphism and reuse). Libmapper has similar properties, but its distributed nature also natively supports collaborative design, since an arbitrary number of GUIs can interact with the same mapping network. It also adds support for supervised machine learning tools through the ability to query the value of any signal in the network, including "destination" signals (synth or application inputs).

We have built a prototype bridge between the two tools in order to exploit their complementarity: special ICon devices can give access to any libmapper signal available on the network. This enables the design of "compound devices" that include signals from a variety of sources gathered together into a logical collection (Fig. 2). By leveraging the benefits of the two approaches, especially their support for the interaction models discussed above, we envision the design of musical applications that treat mappings and interaction as first-class objects.

Paper Substrates and PaperComposer

PaperComposer [START_REF] Garcia | PaperComposer: Creating Interactive Paper Interfaces for Music Composition[END_REF] is an "interface builder" enabling composers to create, manage and use their own interactive paper interfaces with Paper Substrates.

Paper Substrates [START_REF] Garcia | Interactive Paper Substrates to Support Musical Creation[END_REF] are interactive paper components that support the creation and manipulation of complex musical data through user-defined representations, and a paper-based interface, using components such as note containers, tone networks or graph paper (Fig. 3).

Instead of radically replacing current tools, Paper Substrates extend existing music programming environments that composers are familiar with. Paper Substrates rely on existing data processing and storage functions to provide a new layer of user interaction based on handwritten input. They can be seen as interaction instruments rather than as independent software components: they modify or create application objects by communicating directly with the application.

Communication is based on the OSC protocol over data channels, allowing flexibility and independence. Other substrates can subscribe to these data channels to create chains of components acting on each other's data or characteristics before redirecting the result to a particular application. For example, a composer can Routing parts of the libmapper network through the bridge also allows the use of ICon's library of interaction models and the integration of complex state machines into the mapping configuration.

combine a component for entering pitches that uses symbolic notation with another that defines the pitches' amplitudes with a curve, to create a more complex object.

Composers can use the modular components of PaperComposer to create, customize and use their own paper interfaces for a range of musical creation tasks, such as drawing control curves or writing musical sequences with music notation (Fig. 4). New paper components can be developed and integrated within PaperComposer using a Java API, and then used and re-used by composers with any compatible application, supporting co-adaptation.

Conclusion and Future Work

The tools introduced above are examples of technology that embeds knowledge. By bringing together models and tools from HCI and computer music, they facilitate the exploration and definition of new interactive technologies for musical creation and performance. They have been used by performers, composers and instrument designers, and to support the creation of public performances. We are now working on the next generation of our design workbench, which will include more concrete and practical design guidelines. For example, libmapper could suggest signal connections based on well-known HCI theories and models such as integrality and separability of input devices or bimanual interaction. ICon could provide visual feedback on the properties and "quality" of mappings based on similar theories and models, and support more advanced visual programming tools, e.g. for specifying how a Paper Substrate should interpret pen input.

Figure 1 :

 1 Figure 1: Rasmussen's model applied to musical interaction: the three performance behaviors are: skill-based, rule-based or knowledge-(or model)-based (left column). Performance contexts are shown on the right column. At the bottom, contexts demand close temporal coupling between performer and instrument, with little tolerance for interruption. At the top, contexts have much looser coupling.

Figure 2 :

 2 Figure 2: A compound device is created in ICon from various libmapper signals coming from different physical devices. This enables the construction of instrument models from highlevel components; here a device consists of two mobile phones and a Kinect depth camera.Routing parts of the libmapper network through the bridge also allows the use of ICon's library of interaction models and the integration of complex state machines into the mapping configuration.

Figure 3 :

 3 Figure 3: PaperComposer interface. (1) Current document pages. (2) Toolbar to manage document and set parameters. (3) Thumbnails of available component classes can be dragged into the virtual page to create new instances. (4) Virtual page with components: musical staves, play-bar, button, knob slider, tonnetz, curve container.

Figure 4 :

 4 Figure 4: Interacting on paper to explore computer-based musical processes. Photo H. Raguet © Inria.