Fabio Morbidi

Alain Y Kibangou

A Distributed Solution to the Network Reconstruction Problem

Keywords: Network reconstruction, Consensus algorithms, Distributed estimation, Networked systems

 that by collecting noise-contaminated time series generated by a coupled-oscillator system at each node of a network, it is possible to robustly reconstruct its topology, i.e. determine the graph Laplacian. Restricting ourselves to linear consensus dynamics over undirected communication networks, in this paper we introduce a new dynamic average consensus least-squares algorithm to locally estimate these time series at each node, thus making the reconstruction process fully distributed and more easily applicable in the real world. We also propose a novel efficient method for separating the off-diagonal entries of the reconstructed Laplacian, and examine several concepts related to the trace of the dynamic correlation matrix of the coupled single integrators, which is a distinctive element of our network reconstruction method. The theory is illustrated with examples from computer, power and transportation systems.

Introduction

In graph theory the "basic inverse problem" consists in determining the topology of a graph from its adjacency spectrum [START_REF] Van Dam | Developments on spectral characterizations of graphs[END_REF]. It is well known that: i) graphs with few distinct eigenvalues tend to have some kind of regularity; ii) cospectral graphs have the same number of closed walks of a given length; iii) there exist cospectral graphs which are not isomorphic. Therefore, the adjacency (and analogously the Laplacian) spectrum of a graph is not sufficient, in general, to unambiguously identify the graph topology. Recently, the interest in the basic inverse problem has been revitalized by several works in the control and mathematical physics literatures, dealing with the reconstruction of the topology of a network of dynamical systems (this is also sometimes referred to as network "identification" or "exploration" problem). Uncovering the relationship between dynamics and network structure has indeed relevant applications in biology (biochemical, neural and ecological networks), finance, computer science (Internet and World Wide Web), transportation (delivery and distribution networks), and electrical engineering (power grids). In [START_REF] Nabi-Abdolyousefi | Network identification via node knockout[END_REF], the interaction geometry among a known number of agents adopting (weighted) consensus-type algorithms for their coordination, is reconstructed using a grounding procedure inspired by experimental biology, ⋆ Corresponding author, Tel. : +33 (0)4 76 61 53 22 Email addresses: fabio.morbidi@inria.fr (Fabio Morbidi ⋆ ), alain.kibangou@ujf-grenoble.fr (Alain Y. Kibangou).

called "node knockout" (see also [START_REF] Suzuki | Knock-Out/Down-Based Structure Identification in Networks of Heterogeneous Subsystems[END_REF]), while in [START_REF] Kibangou | Decentralized Laplacian Eigenvalues Estimation and Collaborative Network Topology Identification[END_REF] the consensus matrix is determined from its eigenstructure estimated in a distributed fashion. Stochastic methods have lately emerged as powerful alternatives to the deterministic ones in [START_REF] Nabi-Abdolyousefi | Network identification via node knockout[END_REF][START_REF] Suzuki | Knock-Out/Down-Based Structure Identification in Networks of Heterogeneous Subsystems[END_REF][START_REF] Kibangou | Decentralized Laplacian Eigenvalues Estimation and Collaborative Network Topology Identification[END_REF]: the main idea here is to reconstruct the network topology from noise-contaminated time series collected at each node of the graph. Noise is ubiquitous, e.g., in biological networks and relying on such a natural variability as a noninvasive network-identification tool appears extremely promising [START_REF] Hayden | Network Reconstruction from Intrinsic Noise[END_REF]. Several stochastic reconstruction methods have been proposed in the recent literature. In [START_REF] Shahrampour | Topology Identification of Directed Dynamical Networks via Power Spectral Analysis[END_REF], the topology of a directed weighted network of LTI systems is estimated via power spectral analysis, while in [START_REF] M A T E R A S S I | Model identification of a network as compressing sensing[END_REF] the authors formulate the problem of network reconstruction as a compressing sensing problem. Other approaches have exploited the sparsity of the network and used the Bayesian information criterion to measure the graph structure from stationary time series [START_REF] Bach | Learning graphical models for stationary time series[END_REF] or optimizationbased methods [START_REF] Julius | Genetic network identification using convex programming[END_REF]: however, although the sparsity assumption is well justified in some applications (e.g. in biological networks), it may lead to poor results in other cases as shown in [START_REF] Gonçalves | Necessary and sufficient conditions for dynamical structure reconstruction of LTI networks[END_REF]. More related to the standard system-identification literature is the work in [START_REF] Hayden | Network Reconstruction from Intrinsic Noise[END_REF], where the authors formulated the network reconstruction problem as a variant of the spectral factorization problem, and [START_REF] Dankers | Dynamic network identification using the direct prediction-error method[END_REF] where the classical "direct method" of closedloop identification is utilized. However, in these works unknown noise sources are applied only to the states that are measured, which is an unrealistic assumption in many applications. Finally, in [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF], an interesting connection between dynamic correlation and topology in noisy coupled-oscillator networks (where the noise magnitude is known) has been unveiled. The reconstruction algorithm originally proposed in [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF] has been recently improved in [START_REF] Fioriti | Exact Network Reconstruction from Consensus Signals and One Eigenvalue[END_REF] by defining a different thresholding mechanism based on the largest eigenvalue of the Laplacian matrix estimated in a distributed fashion. However, all the aforementioned stochastic approaches are centralized or semi-centralized (as [START_REF] Fioriti | Exact Network Reconstruction from Consensus Signals and One Eigenvalue[END_REF]), since they require the knowledge of time series at each node of the network, and thus impractical for real-world applications. This paper builds upon [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF] and extends it in several new directions. In particular, the original contributions of this work are threefold:

• The approach in [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF] is made fully distributed by leveraging a new dynamic average consensus least-squares algorithm for the local estimation of the noisy time series at each node of the graph. In this way, each node is able to infer the topology of the overall undirected network (cf. Fig. 1): to the best of our knowledge, together with [START_REF] Kibangou | Decentralized Laplacian Eigenvalues Estimation and Collaborative Network Topology Identification[END_REF], this is the only completely-decentralized reconstruction algorithm available in the literature for networks of dynamical systems. An original stability analysis of the proposed dynamic average consensus least-squares estimator is also performed. • A new robust and effective mechanism for separating the -1 and 0 off-diagonal entries of the reconstructed graph Laplacian is proposed. The method is based on κ-means clustering [START_REF] Lloyd | Least-squares quantization in PCM[END_REF] and it is simpler than the thresholding criteria utilized in [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF][START_REF] Fioriti | Exact Network Reconstruction from Consensus Signals and One Eigenvalue[END_REF]. • Interesting connections are shown between the trace of the dynamic correlation matrix of the coupled singleintegrator nodes (which plays a crucial role in our reconstruction mechanism), the total effective resistance of the network [START_REF] Ghosh | Minimizing effective resistance of a graph[END_REF], and the H 2 norm of noisy reduced consensus dynamics [START_REF] Xiao | Distributed average consensus with least-mean-square deviation[END_REF][START_REF] Young | Robustness of noisy consensus dynamics with directed communication[END_REF].

Note that the distributed network-reconstruction algorithm proposed in this paper may represent a valid alternative (at least at small scales) to "web crawlers" for the World Wide Web and to "traceroute sampling" for the Internet, where there does not exist a vantage point with complete information about the overall structure of the system. Moreover, it may be useful for probing the structure of dynamically-changing networks, e.g. road networks, where links can appear/disappear over time because of accidents or works on the carriageway. The rest of this paper is organized as follows. Sect. 2 presents some background material. The main theoretical results of the work are introduced in Sect. 3 and Sect. 4, and the theory is illustrated with numerical simulations on realistic networks in Sect. 5. Finally, Sect. 6 summarizes the main contributions of the paper and outlines some promising future research directions.

Preliminaries

In this section we recall some notions of algebraic graph theory and robust control, and introduce the notation. Let G =( V, E) be an undirected graph (or network), where V = {1,...,n} is the set of nodes and E the set of links. N (i) will indicate the set of nodes adjacent to node i in the graph G. All graphs in this paper are finite, connected, with no self-loops and multiple links.

The adjacency matrix A =[ a ij ]o fg r a p hG is an n × n matrix defined as

a ij = 1i f{i, j}∈E,
0o t h e r w i s e .

The Laplacian matrix L =[ ℓ ij ]o fg r a p hG is an n × n symmetric positive semidefinite matrix defined as L = D -A where D =d i a g ( A1) is the degree matrix and 1 is a column vector of n ones. From this definition, it follows that:

ℓ ij = ⎧ ⎪ ⎨ ⎪ ⎩ n j=1 a ij if i = j, -1i f {i, j}∈E,i , j ∈{1,...,n},
0o t h e r w i s e .

Let L † be the (Moore-Penrose) pseudoinverse of the Laplacian matrix, J = 1 n 11 T the n × n (rank one) averaging matrix, and I n the n × n identity matrix. Then, we have that

L † L = I n -J, (1)
which is the projection matrix onto the image of L [START_REF] Ghosh | Minimizing effective resistance of a graph[END_REF]. From (1), it can be verified that:

L † =( L + J) -1 -J, (2)
and that trace(

L † )= n i=2 1 λ i (L) , (3)
where 0 = λ 1 (L) <λ 2 (L) ≤ ... ≤ λ n (L) are the ordered eigenvalues of L [START_REF] Ghosh | Minimizing effective resistance of a graph[END_REF]. Note that L † inherits from L the property of being symmetric and positive semidefinite. Moreover, L † and L share the same null space. The H 2 norm of a general LTI system ẋ = Ax+ Bu, y = Ex, with A Hurwitz, is given by (trace(EXE T )) 1/2 where the positive semidefinite matrix X (the controllability Gramian) solves the algebraic Lyapunov equation AX + XA T + BB T = 0. In addition, if the pair (A, B) is controllable, X is positive definite. Notation: |S | will denote the cardinality of the set S , ∅ the empty set, Re(z)a n d|z| the real part and modulus of the complex number z, respectively, ⊗ the Kronecker product, blkdiag(•) a block-diagonal matrix, E[ • ] the expectation operator, • 2 the Euclidean norm of a vector and • L2 the total energy or L 2 norm of a vector-valued signal.

Problem formulation

For the reader's convenience we briefly review here the main results in [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF], which form the basis for our subsequent developments. We will start with a general formulation dealing with directed graphs, and then we will specialize our results to undirected networks with coupled single-integrator node dynamics. Consider a directed network of n nonidentical coupled oscillators where x i ∈ IR m , i ∈{ 1,...,n}, denotes the state variable of the i-th oscillator and x i (0) the initial state. In the presence of noise, the dynamics of the whole oscillator system can be expressed as

ẋi = F i (x i ) -γ n j=1 ℓ ij H(x j )+η i , (4)
where

F i :I R m → IR m is the intrinsic dynamics of the i-th oscillator, H :I R m → IR m
is the coupling function of the oscillators, γ>0 is the coupling strength, and η i ∈ IR m is the zero-mean white noise with covariance matrix Q i associated with the i-th oscillator. We will assume that x i (0) and η i are uncorrelated for any i ∈ {1,...,n},an dt h atη i and η j are uncorrelated for any i = j, i, j ∈{ 1,...,n}. Note that (4) encompasses several well-known coupled-oscillator systems, such as Kuramoto phase oscillators and Rössler dynamics. Let x i be the counterpart of x i in the absence of noise and ξ i a small perturbation. Then we can write x i = x i + ξ i and plugging it in (4), we obtain the variational equation

ξi = J[F i (x i )] ξ i -γ n j=1 ℓ ij J[H(x j )] ξ j + η i ,
where J[F i (•)] and J[H(•)] are the m × m Jacobian matrices of F i and H, respectively. Now, let ξ =[ ξ T 1 ,...,ξ T n ] T ∈ IR mn be the deviation vector, η =[η T 1 ,...,η T n ] T ∈ IR mn the noise vector, and

C 1 T T 0 ξ(t) ξ T (t) dt, (5)
the dynamic correlation matrix of the oscillators, where [0, T] is the (long) time interval in which the time average is performed. We then have that:

0 = d dt 1 T T 0 ξ(t) ξ T (t) dt = MC+ CM T + 1 T T 0 η(t) ξ T (t) dt + 1 T T 0 ξ(t) η T (t) dt, (6)
with

M = J[F(x)] -γ L ⊗ J[H(x)] where J[F(x)] = blkdiag(J[F 1 (x 1 )],...,J[F n (x n )]) and x =[x T 1 ,...,x T n ] T ∈ IR mn .
Assuming that M is Hurwitz, after few algebraic manipulations we can rewrite [START_REF] Hayden | Network Reconstruction from Intrinsic Noise[END_REF] as

MC + CM T +2Q = 0, (7)
where

Q = blkdiag(Q 1 ,...,Q n ).
Note that [START_REF] Shahrampour | Topology Identification of Directed Dynamical Networks via Power Spectral Analysis[END_REF], an algebraic Lyapunov equation in the unknown C, reveals a general relationship between the dynamic correlation matrix C and the Laplacian L in the presence of noise with covariance matrix Q.Itisw ellkno wnthatif M is Hurwitz and Q positive semidefinite, the unique solution of ( 7) is positive semidefinite and it has the representation:

C =2 ∞ 0 e M t Q e M T t dt. (8)
If we assume that M = M T and Q = σ 2 2 I mn , σ>0, the integral in (8) can be explicitly computed and we obtain C = σ 2 2 M -1 . If we further assume a linear coupling between the oscillators, i.e. J[H(•)] = I mn ,n oi n t r i n s i c dynamics, i.e. F i (•)=0 for all i, and one-dimensional state variables, i.e. m = 1, we find the following simple formula relating the dynamic correlation matrix C to the Laplacian L:

C = σ 2 2 γ L † . (9)
Note that this corresponds to the case of nodes with coupled single-integrator dynamics driven by white noise, communicating over an undirected graph G, i.e. to the noisy consensus dynamics

ẋ = -γ Lx + η, (10)
where x =[ x 1 ,...,x n ] T ∈ IR n and η ∈ IR n . This simplified yet meaningful scenario will be considered in the rest of this paper. Henceforth, we also assume that the number n of nodes of G is known. Note that from (9), we obtain

L = σ 2 2 γ C † = 2 γ σ 2 C + J -1 -J, (11)
where the second equality follows from (2). We then see that the entire network structure, embodied in the Laplacian matrix L, can be extracted from the knowledge of noise-corrupted time series. In fact, recalling the definition of ξ, the dynamic correlation matrix can be computed as

C = 1 T T 0 (x(t) -x(t)1)(x(t) -x(t)1) T dt,
where

x(t)= 1 n n i=1 x i (t)= 1 n 1 T x(t).
Since x(t)x(t)1 =(I n -J) x(t), we can rewrite equation [START_REF] Gonçalves | Necessary and sufficient conditions for dynamical structure reconstruction of LTI networks[END_REF] as

L = 2 γ σ 2 (I n -J) C R (I n -J)+J -1 -J, (12)
where

C R = 1 T T 0 x(t) x T (t) dt. ( 13
)
Since the integral in ( 13) cannot be analytically computed in practice, assuming S discrete observations in t h et i m ei n t e rv al[ 0, T], we can approximate C R as

C R ≃ 1 S S-1 k =0 x(k) x T (k). ( 14
)
Note that in order to compute the right-hand side of equation ( 14), we need to know the time series {x i (0),x i (1), ..., x i (S -1)} at each node of the graph. This means that in its present form, the network reconstruction algorithm relies on global information and it is thus centralized.

Next, we will introduce a practical method to separate the off-diagonal entries of the reconstructed Laplacian (Sect. 3.1), and highlight parallels between the trace of the dynamic correlation matrix C and concepts from other fields (Sect. 3.2). In Sect. 4, we will finally present a consensus-based least-squares estimation algorithm in order to make the reconstruction process fully distributed and computationally treatable:i nt h i sw a y ,b y only interacting with its neighbors, each node will be able to locally and efficiently estimate the topology of the overall network.

Separation of the off-diagonal entries of the reconstructed Laplacian

Fig. 2 shows the distribution of the off-diagonal entries of matrix σ 2 2γ C † for a social network with 62 nodes and 159 links [START_REF] Lusseau | The bottlenose dolphin community of Doubtful Sound features a large proportion of long-lasting associations[END_REF], obtained by simulating system [START_REF] Julius | Genetic network identification using convex programming[END_REF] with γ =1a n dσ 2 = 2 for 2000 seconds from a random initial condition. The distribution has one peak centered at -1, corresponding to the existent links in the network and the other centered at 0, corresponding to the zero elements of L.N o t et h a ti fL were perfectly reconstructed from σ 2 2γ C † , the two peaks would be sharp vertical lines. Since this does not occur in practice, a threshold needs to be determined to separate the two peaks 2γ C † for a network with 62 nodes [START_REF] Lusseau | The bottlenose dolphin community of Doubtful Sound features a large proportion of long-lasting associations[END_REF]. We set γ =1 ,σ 2 =2a nd run [START_REF] Julius | Genetic network identification using convex programming[END_REF] for 2000 seconds. [Assignment step (Forgy )]: Assign each c h to the cluster S κ (q) whose center µ κ (q) is closest (in the squared Euclidean-distance sense) to c h .

4:

[Update step]: For κ ∈{1, 2}, update the centers µ κ (q)as

µ κ (q +1) = 1 |S κ (q)| c h ∈ Sκ(q) c h .
5: end for and distinguish the existent from the nonexistent links. An analytical procedure was proposed in [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF] for finding such a threshold: however, it is neither immediate, nor appears well-physically motivated. On the other hand, the thresholding mechanism in [START_REF] Fioriti | Exact Network Reconstruction from Consensus Signals and One Eigenvalue[END_REF], which is based on the largest eigenvalue of the graph Laplacian, is not described in sufficient detail sot h a ti tc a nb er e p l i c a t e d herein.

As a simple yet effective alternative to the thresholding methods in [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF][START_REF] Fioriti | Exact Network Reconstruction from Consensus Signals and One Eigenvalue[END_REF], in this paper we propose to use κ-means clustering [START_REF] Lloyd | Least-squares quantization in PCM[END_REF]. Let c =[ c h ] ∈ IR n(n-1) be the vector of off-diagonal entries of σ 2 2γ C † listed in the natural order. In our case, the clustering problem amounts to partition the entries of c into two clusters S 1 , S 2 where S 1 is relative to the -1 values and S 2 to the 0 values (S 1 ∩ S 2 = ∅, S 1 ∪ S 2 = {c 1 ,c 2 ,...,c n(n-1) }), so as to minimize the squared distortion:

min S1, S2,µ1,µ2 κ ∈{1,2} c h ∈ Sκ 1 2 c h -µ κ 2 ,
where µ κ is the center of cluster S κ , κ ∈{ 1, 2}.T h i s combinatorial optimization problem is nonconvex and exhaustive search for the optimal clusters is intractable. This intractability motivates low-complexity heuristics such as the κ-means algorithm. The algorithm uses an iterative refinement technique: in fact, it alternates between minimizing the squared distortion over {S 1 , S 2 } (assignment step) and {µ 1 ,µ 2 } (update step), see Algorithm 1. The algorithm has converged when the assignments no longer change. Since both the steps optimize the squared distortion and there only exists a finite number of assignments, the algorithm must converge to a (local) minimum after N < ∞ steps. Once convergence has been achieved, the Laplacian matrix can be recovered, for h ∈{ 1,...,n(n -1)}, i = j, i, j ∈{1,...,n},as

ℓ ij = ⎧ ⎨ ⎩ -1i fc h ∈ S 1 (N -1), 0i f c h ∈ S 2 (N -1),
and

ℓ ii = n j=1,j =i | ℓ ij |, i ∈{1,.
..,n}.

Physical interpretations of the correlation matrix

Effective graph resistance

Assume the graph G is seen as an electrical circuit, where a link corresponds to a unit resistor. The effective resistance between a pair of nodes i and j, is the electrical resistance measured across nodes i and j.T h etotal effective resistance, R tot , is defined as the sum of the effective resistance between all distinct pairs of nodes of G [START_REF] Ghosh | Minimizing effective resistance of a graph[END_REF]. It has been shown (see, e.g., [15, Sect. 2.5]) that R tot depends on the spectrum of the graph Laplacian as

R tot = n n i=2 1 λ i (L)
.

Hence, by recalling (3), we obtain the following concise formula relating the total effective resistance and the dynamic correlation matrix C R tot = 2 γn σ 2 trace(C).

Consensus robustness and the H 2 norm

Since system ( 10) is only marginally stable in the noisefree case, let us restrict to the dynamics on span(1) ⊥ , i.e. on the subspace of IR n orthogonal to the subspace spanned by 1.L e tt h em a t r i xΘ ∈ IR (n-1) × n be such that Θ 1 = 0, ΘΘ T = I n-1 and define r Θx. By differentiating r with respect to time, we obtain:

ṙ = -γ L * r + Θ η, (15)
where L * = ΘLΘ T is the reduced Laplacian matrix (note that L * has the same eigenvalues of L except for the zero eigenvalue). Let us assume that the output equation of system ( 15) is y = I n-1 r, (16) then we get the following result (for the proof the reader is referred to [17, Lemma 2]).

Proposition 1 Let the undirected graph G be connected and σ 2 =1.TheH 2 norm of system ( 15)-( 16) is

H = n i=2 1 2 γλ i (L) 1/2 = (trace(C)) 1/2 .

Let us now define Σ(t)

E[r(t) r T (t)]. Note that H in Prop. 1 coincides with the measure of robustness to white-noise inputs of system (15) considered in [START_REF] Young | Robustness of noisy consensus dynamics with directed communication[END_REF] H

lim t→∞ E[ r(t) 2 ] = lim t→∞ (trace(Σ(t))) 1/2 .
This definition of robustness corresponds to the steadystate least-mean-square deviation introduced in [START_REF] Xiao | Distributed average consensus with least-mean-square deviation[END_REF] for discrete-time consensus dynamics: a related H 2 normbased notion quantifying the deviation from consensus, is called (first-order) network coherence in [START_REF] Bamieh | Coherence in large-scale networks: Dimension-dependent limitations of local feedback[END_REF]. Finally, under the same assumptions of Prop. 1, we have the following geometric interpretation of matrix C in relation to the properties of the controllability Gramian (2 γ L * ) -1 of system [START_REF] Ghosh | Minimizing effective resistance of a graph[END_REF], cf. [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF]Th. 6.18].

Proposition 2 (Reachable set with bounded noise)

The closure of the set of all the states of system (15) reachable from r(0) = 0 with noise bounded as

η 2 L2 ∞ 0 η T (t) η(t) dt ≤ 1,
is given by the ellipsoid

T = r : p T r ≤ (p T ΘCΘ T p) 1/2 , ∀ p ∈ IR n-1 .
In other words, Prop. 2 says that the points in the state space of system (15) that can be reached using an excitation with total energy one, are given by an ellipsoid determined by C.

Distributed network-state estimation

As aforementioned, in this section we will present a consensus-based strategy in order to make the reconstruction procedure described in Sect. 3 fully distributed.

To this end, let us introduce the following measurement equation for node i in the network G:

y i (k)=E i x(k),k ∈{0, 1,...,S -1},
where E i ∈{ 0, 1} (|N (i)|+1) × n is a selection matrix (the number of observations S is assumed to be known a priori by each node). In other words, at time instant k ∈{0, 1,...,S-1},nodei measures its own state x i (k) (the first component of vector y i (k)), and the state x j (k)ofj ∈N(i). By collecting the measurements of all the agents at time k,w eo b t a i ny(k)=Ex(k), k ∈{ 0, 1,...,S -1},w h e r ey =[ y T 1 ,...,y T n ] T and E =[ E T 1 ,...,E T n ] T . If the columns of E are linearly independent (i.e. E T E is invertible), then the leastsquares estimate of x(k)i s

x(k)=E † y(k)=( E T E) -1 E T y(k) = 1 n n i=1 E T i E i -1 1 n n i=1 E T i y i (k) . ( 17
)
Note that the two quantities inside the square brackets in [START_REF] Young | Robustness of noisy consensus dynamics with directed communication[END_REF], can be computed in a distributed fashion using consensus-based algorithms. Let us then assume that node i maintains a matrix P i ∈ IR n×n and a vector z i ∈ IR n , and executes the following iterations for h ∈ {0, 1, ..., D-1}, D>1,

P i (h +1) = P i (h)+ς j ∈N(i) (P j (h) -P i (h)), ( 18a
)
z i (h +1) = z i (h)+ς j ∈N(i) (z j (h) -z i (h)), ( 18b
)
which are respectively initialized with

P i (0) = E T i E i , z i (0) = E T i y i (k), i ∈{ 1,.
..,n}. If the step size ς ∈ (0, 1) satisfies the condition 1 ς<2/λ n (L), cf. [START_REF] Mesbahi | Graph Theoretic Methods in Multiagent Networks[END_REF]Lemma 8.1], then lim

h→∞ P i (h)= 1 n n i=1 E T i E i , lim h→∞ z i (h)= 1 n n i=1 E T i y i (k)
, and node i asymptotically obtains an estimate x (i) of vector x at time k according to:

x (i) (k) = lim h→∞ P -1 i (h) z i (h). (19)
Node i can then compute its local estimate

C (i) R of C R as C (i) R = 1 S S-1 k =0 x (i) (k)( x (i) (k)) T ,
and eventually obtain an estimate of the Laplacian matrix via equation ( 12) and the application of κ-means clustering. Note that the iterates P i (h)m a ynot be invertible for all values of h: hence, estimate [START_REF] Bamieh | Coherence in large-scale networks: Dimension-dependent limitations of local feedback[END_REF] at node i can only be computed once P i (h) becomes nonsingular.

We also observe here that differently from the iterations in (18a), which have to be performed only once, the iterations in (18b) must be executed at each time instant k ∈{0, 1, ...,S -1}, which is computationally demanding, especially for large networks. In order to reduce the computational load, we can replace (18b) with a dynamic average consensus estimator which tracks the average of all the nodes' inputs E T i y i (k)f o r k ∈{0, 1,..., S -1}. Because of its attractive robustness properties, we will make use of the discrete-time version of the PI dynamic average consensus estimator introduced in [START_REF] Lynch | Decentralized environmental modeling by mobile sensor networks[END_REF]:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ z i (k +1)=(1-εα) z i (k)+ε K P j ∈N(i) (z j (k) -z i (k)) -ε K I j ∈N(i) (w j (k) -w i (k)) + εαE T i y i (k), w i (k +1)=w i (k)+ε K I j ∈N(i) (z j (k) -z i (k)), (20)
where z i is node i's estimate of the average of all the nodes' inputs, w i is an internal estimator state, K P , K I > 0 are estimator gains, ε>0 is the step size, and α>0 is a parameter governing the rate at which new information replaces old information in the dynamic averaging process (the so-called "forgetting factor"). Each node communicates its estimate z i and its internalestimator state w i to its neighbors in the network. Note that differently from a static average consensus protocol, in [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF] the inputs continually drive the estimator and each z i (k) approximately tracks the true average of

E T i y i (k). By defining z =[ z T 1 ,...,z T n ] T ∈ IR n 2 and w =[ w T 1 ,...,w T n ] T ∈ IR n 2
, we can rewrite system (20) in compact form as:

z(k +1) w(k +1) =(Γ ⊗ I n ) z(k) w(k) + εαI n 2 0 E T y(k), (21)
the overall network by each node. However, it is well-known that λn(L) ≤ n [START_REF] Godsil | Algebraic graph theory[END_REF], and the upper bound n c a nbeu s e db y each node to determine ς in [START_REF] Lusseau | The bottlenose dolphin community of Doubtful Sound features a large proportion of long-lasting associations[END_REF].

where

Γ = (1 -εα)I n -ε K p L ε K I L -ε K I LI n .
The proof of the next proposition, which is not available in the literature, is given in the Appendix.

Proposition 3 (Stability of the PI estimator)

The PI estimator ( 21) is stable if the step size

ε< 2Re(β i ) |β i | 2 , ∀ i ∈{2,...,n}, (22)
where

β i = 1 2 α -K P ζ i ± (α -K P ζ i ) 2 -4K 2 I ζ 2 i , (23)
and ζ i , i ∈{2,...,n},i sth ei-th eigenvalue of -L.

Remark 1 As the size of the network G increases, longer time series are needed to get the same level of reconstruction accuracy (and, in general, to make P i (h) invertible): in fact, the information takes more time to flow through G and more communication rounds are necessary. On the other hand, higher noise magnitude may reduce, in general, the number of necessary consensus steps but also increase the power dissipation at each node. Therefore, computational complexity/energy consumption needs to be weighted against reconstruction accuracy in practical applications where available resources are at premium. ⋄

Numerical validation

Numerical simulations have been conducted with five realistic networks in order to compare the performance of the centralized and distributed network reconstruction methods.

• Data and computer networks: A full binary tree is a rooted tree in which every node other than the leaves has two children [START_REF] Godsil | Algebraic graph theory[END_REF]. Full binary trees are extensively used to organize data in computer systems. In our simulations, we considered the "Full binary tree T 4 "which has 4 levels (excluding the root node) and 31 nodes.

Hypercube networks [START_REF] Godsil | Algebraic graph theory[END_REF] have found frequent applications in parallel processing (e.g. in massively-parallel supercomputers). In this study, we chose the "Hypercube Q 4 ", a regular network with 16 nodes. • Electric power network : the "IEEE 9-bus" power system consists of four generators (nodes from 1 to 4) and five loads (nodes from 5 to 9), see [START_REF] Suzuki | Knock-Out/Down-Based Structure Identification in Networks of Heterogeneous Subsystems[END_REF]. • Transportation networks: the links of the network "Road Grenoble" correspond to the major arterials in the city of Grenoble in France, and the nodes to 27 road intersections (see Fig. 3(a)). The metro of Lyon, "Metro Lyon", in France, has a total of 4 lines and 42 stations [START_REF] Derrible | Network centrality of metro systems[END_REF]. In this study, only the 6 termini (black) and 4 transfer stations (white) are considered as nodes (see Fig. 3(b)).

Following [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF], the network-reconstruction performance was quantified using the success rate of existing and non-existing links (SREL, SRNL, respectively). SREL (SRNL) is defined as the ratio of the number of successfully predicted existent (nonexistent) links to the total number of existent (nonexistent) links. In our tests, the coupling strength γ was set to 1 and the noise magnitude σ 2 = 2. Random initial conditions x(0), and z(0), w(0) were selected for system [START_REF] Julius | Genetic network identification using convex programming[END_REF] and the PI estimator [START_REF] Godsil | Algebraic graph theory[END_REF], respectively. κ-means clustering was used in all cases for separating the off-diagonal entries of the reconstructed Laplacian matrix, and it produced results more consistent and repeatable than those of the thresholding method in [START_REF] Ren | Noise bridges dynamical correlation and topology in coupled oscillator networks[END_REF]. Table 1 reports the SREL/SRNL obtained as the average of 10 trials of 1000 seconds (ς =0.1, S = 10000) for the centralized and distributed algorithms, and the value of the parameters used in the numerical tests (last column). In the distributed case, we reported the SREL/SRNL of two adjacent nodes drawn at random in the network: in most of the cases a success rate greater than 0.9 was obtained. In order to illustrate the performance of our dynamic average consensus least-squares algorithm, Fig. 4 shows the time evolution of x j (red) and x (6) j (black), j ∈{1, 2, 3, 4},obtained via (18a) and ( 21), for the hypercube network Q 4 (note that node 6 is adjacent to node 3 but not to nodes 1, 2, 4). As it is evident from the figure, the estimated state components satisfactorily track the true state.

Conclusions and future work

In this paper we have presented a fully distributed algorithm for reconstructing the topology of an undirected consensus network. Our approach consists in locally estimating the noise-corrupted network state by using a distributed least-squares algorithm. We proposed a new efficient method for separating the off-diagonal entries of the reconstructed graph Laplacian based on κ-means clustering, and offered some insight into the trace of the dynamic correlation matrix of the coupled single integrators. Numerical simulations performed on realistic networks have demonstrated the effectiveness of our distributed reconstruction method. (black), j ∈{1, 2, 3, 4}, for the hypercube network Q4.

In this paper, we have only dealt with the Boolean structure of a network (i.e. with the existence or nonexistence of links between nodes). In future works, we would like to extend the proposed method to recover the weights of a network, and to estimate the directions of the links in directed graphs. We are also going to explore other methods for computing the sample covariance matrix in equation [START_REF] Lloyd | Least-squares quantization in PCM[END_REF] in a distributed fashion (e.g. by estimating its eigenvectors in a decentralized way), and to relax the assumption of known noise magnitude. Finally, it is worth pointing out that because of its computational complexity, the proposed reconstruction algorithm is not applicable, at present, to networks with hundreds or thousands of nodes. In order to address this problem, we plan to exploit the sparsity of large-scale real-world networks for distributed least-squares estimation, and to leverage the information coming from anchor or ancillary nodes around the network.

Appendix

λ ± i =1- ε 2 α -K P ζ i ± (α -K P ζ i ) 2 -4K 2 I ζ 2 i .
Note that λ + 1 =1-εα and λ - 1 =1s i n c eζ 1 =0 ,a n d that λ - 1 is the only eigenvalue of Γ having unitary modulus. Let us now define the complex number β i as in equation ( 23) (notice that Re(β i ) ≥ 0, ∀ i). Then for system [START_REF] Godsil | Algebraic graph theory[END_REF] to be stable, we must impose that |1-εβ i | < 1, ∀ i ∈{2,...,n}, from which condition [START_REF] Mesbahi | Graph Theoretic Methods in Multiagent Networks[END_REF] follows.

Fig. 1 .

 1 Fig. 1. Distributed network reconstruction:n o d ei must reconstruct the topology of the network G by only exchanging information with its set of neighbors N (i)( r e dn o d e s ) .

25Fig. 2 .

 2 Fig. 2. Distribution of the off-diagonal entries of matrix σ 22γ C † for a network with 62 nodes[START_REF] Lusseau | The bottlenose dolphin community of Doubtful Sound features a large proportion of long-lasting associations[END_REF]. We set γ =1 ,σ 2 =2a nd run[START_REF] Julius | Genetic network identification using convex programming[END_REF] for 2000 seconds.

Algorithm 1 :

 1 κ-means clustering 1: [Initialization]: Randomly draw the centers µ κ (0), κ ∈{1, 2},from{c 1 ,c 2 ,...,c n(n-1) }. 2: for q ∈{0, 1,...,N -1} do 3:

Fig. 3 .

 3 Fig. 3. Examples: (a) Network of major arterials in the city of Grenoble (image source: Google maps); (b) Simplified plan of the metro of Lyon (www.tcl.fr).

Fig. 4 .

 4 Fig. 4. Time evolution of xj (red) and x

Lemma 1

 1 Let A ∈ IR n×n and B ∈ IR m×m . Then the eigenvalues of A ⊗ B are of the form λ i (A) λ j (B), i ∈ {1,...,n}, j ∈{1,...,m}. Lemma 2 (Schur's formula) Let A, B, C, D ∈ IR n×n .I fA and C commute, i.e. AC = CA,t h e n det AB CD =d e t ( AD -CB).P r o o fo fP r o p o s i t i o n3 :Owing to Lemma 1, to prove the statement we only need to study the eigenvalues of matrix Γ, which are the roots of det(λI 2n -Γ)=0 .I f we apply Lemma 2, we find thatdet(λI 2n -Γ)= det (λ -1)[(λ -1+εα) I n + ε K P L]+ε 2 K 2 I L 2 . (24) Note now that det(λI n + L)= n i=1 (λζ i )whereζ i is the i-th eigenvalue of -L with 0 = ζ 1 >ζ 2 ≥ ... ≥ ζ n .By comparing[START_REF] Derrible | Network centrality of metro systems[END_REF] with this last equation, we see that det(λI2n -Γ)= n i=1 (λ 2 -(2εα + εK P ζ i )λ +1εα + εK P ζ i + ε 2 K 2 I ζ 2 i) from which it follows that the eigenvalues of Γ can be obtained by solving the equationλ 2 -(2-εα+ε K P ζ i )λ+1-εα+ε K P ζ i +ε 2 K 2 I ζ 2 i =0 . Hence, if λ ±i are the eigenvalues of Γ associated with ζ i , we have that

Table 1

 1 Success rate of existent and non-existent links (SREL/SRNL) with the centralized and distributed reconstruction methods for five different network topologies. The results are the average of 10 trials, each one of 1000 seconds.

	Network	#n o d e s # links SREL/SRNL Centralized SREL/SRNL Distributed (2 nodes)	K P ,K I , α, ε, D
	Full binary tree T 4	31	30	1.0000 1.0000	0.9167 0.9722, 0.8667 0.9733	60, 0.05, 3, 0.005, 2500
	Hypercube Q 4	16	32	1.0000 1.0000	0.9906 0.9977, 0.9127 0.9943	40, 0.05, 8, 0.005, 2000
	IEEE 9-bus	9	9	1.0000 1.0000	0.9000 0.9630, 0.9667 0.9556	50, 0.05, 1.5, 0.005, 2000
	Road Grenoble	27	42	0.9976 0.9997	0.7143 0.8829, 0.7071 0.8903	40, 0.05, 1.5, 0.003, 2000
	Metro Lyon	10	10	1.0000 1.0000	0.7700 0.9971, 0.8000 0.9943	40, 0.05, 1.5, 0.005, 2000

Note that the knowledge of the largest eigenvalue of the Laplacian, λn(L), entails the knowledge of the topology of


Dr Virginie Nazabal
email: virginie.nazabal@univ

P Němec

M Olivier

E Baudet

Optical properties of (GeSe 2 ) 100-x (Sb 2 Se 3 ) x glasses in near-and middle-infrared spectral regions

Keywords: Codes: 42.70.Km, 42.70.Ce, 42.65.-k, 78.20.-e, 81.05.Gc A. glasses, A. chalcogenides, D. optical properties, ellipsometry, mid-infrared

Chalcogenide glasses from Ge-Sb-Se ternary system with different antimony content were fabricated and basic physico-chemical properties (chemical composition, thermal characteristics, density) were evaluated considering the glassy network connectivity. Optical properties of the glasses were heavily studied employing transmission spectroscopy, prism coupling technique, and particularly spectroscopic ellipsometry covering broad spectral range (0.3-20 m). Refractive indices data show very good agreement between ellipsometry and prism coupling techniques in near-infrared range. Moreover, the reliability of infrared spectroscopic ellipsometry was demonstrated for precise determination of refractive index of chalcogenide glasses in near-as well as middle-infrared spectral regions.

Introduction

Chalcogenide glasses, based on S, Se or Te elements in combination with suitable elements from 14. or 15. group of the periodical system (Ge, As, Sb, etc.), are known to present photosensitivity to light exposure [1][2][3][4] and to be highly nonlinear [5,6]. Indeed, they present large Kerr nonlinearities at femtosecond time scale and Their high third order susceptibility χ (3) can reach values several orders of magnitude larger than that of silica glass [7]. For these reasons, chalcogenide glasses have been studied for ultrafast switching in telecommunication systems [8]. Further, chalcogenide glasses can also be used for optical sensing applications from visible to near mid-IR [9][10][11][12]. For mentioned applications, chalcogenide thin films are needed; therefore a key prerequisite is an identification of appropriate bulk glass composition.

The Ge-Sb-Se glassy system is of interest due to its large glass forming region [START_REF] Popescu | Non-Crystalline Chalcogenides[END_REF] giving a possibility of tailoring glass properties in broad ranges. In comparison with chalcogenide glasses compositions based on arsenic, the use of antimony is beneficial regarding the limitation of arsenic in the glass composition in accordance with the recommendations of the current regulation. Moreover, the substitution of As with Sb may increase the (non)linear refractive of the corresponding amorphous chalcogenides due to enlargement of (hyper)polarisability. Further, presence of antimony in amorphous chalcogenides is known to reduce photosensitivity of the material and this property could be tailored by changing the proportion of Sb.

For applications in photonics, one of the most important optical characteristics of chalcogenide glasses is refractive index. In visible and near-infrared (NIR) spectral region, refractive index is often evaluated by modeling of spectroscopic data (transmission and/or reflection), prism coupling method, ellipsometry or it can be calculated by quantum mechanics, etc. [START_REF] Palik | Handbook of Optical Constants of Solids[END_REF]. In this paper, beyond the description of basic physico-chemical properties, we are focused on linear optical properties of selected bulk chalcogenide glasses from Ge-Sb-Se system. This is a required step in order to perfectly master the optical properties of the pseudo-binary (GeSe 2 ) 100-x (Sb 2 Se 3 ) x system for applications whether in the field of optical sensors or in non-linear optics. Special attention was paid to the use of spectroscopic ellipsometry in broad spectral range with the aim to determine linear refractive index from ultraviolet to middle-infrared (MIR).

Material and methods

Bulk samples from pseudo-binary (GeSe 2 ) 100-x (Sb 2 Se 3 ) x system (5  x  70) were synthesized from elements (Ge, Sb, and Se) of high purity (5N) using the conventional melting and quenching technique. After the quenching, glass rods were annealed (20 °C below their glass transition temperature (T g ), 6 hours) and slowly cooled down to room temperature. For optical characterization, resulting glass rods (25 mm diameter) were sliced and polished. A scanning electron microscope with an energydispersive X-ray analyzer (EDS, JSM 6400-OXFORD Link INCA) was used for the determination of Ge-Sb-Se glasses chemical composition and its uniformity. Transmission spectra were recorded in VIS-NIR spectral range with 1 mm thick samples (Perkin Elmer Lambda 1050 spectrophotometer). After calculation of absorption coefficients, cut-off wavelengths were determined for an absorption coefficient of 10 cm -1 . Linear refractive indices of fabricated glasses were measured using the prism coupling technique (Metricon 2010) at 1311 and 1551 nm. Silicon prism was used for these measurements. Linear refractive indices spectral dependencies as well as optical band gap values of (GeSe 2 ) 100-x (Sb 2 Se 3 ) x glasses were obtained from the analysis of variable angle spectroscopic ellipsometry (VASE) data measured using an ellipsometer with automatic rotating analyzer (J. A. Woollam Co., Inc.). The measurement parameters are as follows: spectral region 300-2300 nm (UV-NIR) with wavelengths steps of 20 nm, angles of incidence 50°, 60°, and 70°. For the analysis of VASE data, we used Cody-Lorentz model which includes both the correct band edge function and weak Urbach absorption tail [START_REF] Cody | Semiconductors and semimetals[END_REF]; this model is appropriate for the description of amorphous chalcogenides optical functions [4,[START_REF] Nemec | [END_REF]. Linear refractive indices were extracted also from spectroscopic ellipsometry data measured in NIR-MIR spectral region with IR-VASE ellipsometer (J. A. Woollam Co., Inc.) exploiting rotating compensator. NIR-MIR ellipsometric data were measured within 1700-20000 nm range using wavelength steps and angles of incidence as above. To derive refractive indices in NIR-MIR spectral range, Cauchy like dispersion relation was exploited, setting extinction coefficient in first approximation to zero.

Results and discussion

The bulk (GeSe 2 ) 100-x (Sb 2 Se 3 ) x samples were amorphous as confirmed by X-ray diffraction data, excluding (GeSe 2 ) 30 (Sb 2 Se 3 ) 70 material, which was (at least partly) crystalline and was not used for further optical characterization. As determined using EDS, chemical composition of fabricated glasses is in good agreement with nominal one; the differences are about 1 at. % which correspond to the EDS measurements uncertainty.

No crystallization peaks were observed in DSC curves except for glass with x=70 which is close to glassy domain edge. T g decreases monotonously when antimony concentration increases (from ~360 °C for (GeSe 2 ) 95 (Sb 2 Se 3 ) 5 glass to ~230 °C for (GeSe 2 ) 40 (Sb 2 Se 3 ) 60 composition). This trend follows changes of vitreous network connectivity occurring when antimony is incorporated in chalcogenide glasses. We note that antimony coordination is 3 (supposing [SbSe 3/2 ] pyramids as basic structural motifs containing antimony atoms), whereas germanium is known to be in a tetrahedral environment ([GeSe 4/2 ] tetrahedra) [17][18][19]. On the other hand, density of (GeSe 2 ) 100-x (Sb 2 Se 3 ) x glasses clearly increases with rising x (from ~4.35 g.cm -3 for (GeSe 2 ) 95 (Sb 2 Se 3 ) 5 to ~5.24 g.cm -3 for (GeSe 2 ) 40 (Sb 2 Se 3 ) 60 ). Considering atomic masses of elements (A r (Ge) = 72.64, A r (Sb) = 121.76, A r (Se) = 78.96), the density of glasses increases coherently with Sb 2 Se 3 introduction [20]. This also tends to prove that the substitution of GeSe 2 by Sb 2 Se 3 (leading to a lower degree of connectivity) does not increase drastically the free volume of the glassy network which should have an antagonistic effect on the density.

Optical transmission of studied glasses reaches ~64% in the NIR region; this limit comes from Fresnel reflections at air/glass interface. Fig. 1 gives spectral dependences of absorption coefficient calculated from measured transmission spectra. One can see the shift of the fundamental absorption edge to lower energies with increasing content of antimony in the glasses. The influence of atomic orbitals originated from antimony selenide in the band structure of the selenide glass can be easily felt on the optical band-gap. Cut-off wavelengths (described above) of all the glassy samples are given in Table 1 showing similar tendency.

Linear refractive indices of (GeSe 2 ) 100-x (Sb 2 Se 3 ) x glasses measured at 1311 and 1551 nm by the prism coupling technique are listed in Table 1. As anticipated, the refractive index rises particularly by increasing the proportion of antimony (from 2.41 for (GeSe 2 ) 95 (Sb 2 Se 3 ) 5 glass to 3.01 for (GeSe 2 ) 40 (Sb 2 Se 3 ) 60 composition, at 1550 nm). Consequently, one can expect a similar trend for nonlinear refractive index [21]. Spectroscopic ellipsometry data analysis with Cody-Lorentz model provides a determination of refractive indices dispersion over broad spectral range (300-2300 nm) including absorption onset and fundamental absorption of studied chalcogenide glasses, with an uncertainty of ±1.10 -2 . The Cody-Lorentz model refractive indices results, illustrated in Fig. 2 and Table 1, are in pretty good agreement with the results obtained by the prism coupling technique even thought that absolute precision is higher for prism coupling (±1.10 -3 ). Optical band gap values, extracted from VASE data again employing Cody-Lorentz model, are listed in Table 1 confirming the compositional tendency of cut-off wavelengths even if absolute values of cut-off energies differ from VASE band gap energies due to different method of calculation.

To date, published refractive index data of chalcogenide glasses/thin films in MIR spectral region are infrequent [22]. This can easily be explained by the fact that apart from traditional glasses, these measurements require specific manufacturing of chalcogenide prism involving perfect homogeneity and appropriate size as well as access to an accurate optical metrology set-up. Also, the fabrication of optical components for MIR is often empirical because their dispersion of refractive indices in MIR have not necessarily been characterized enough accurately. Our data show applicability of spectroscopic ellipsometry for the determination of refractive indices of chalcogenide glasses covering 1700-20000 nm region with simple model used (Fig. 2) for a series of (GeSe 2 ) 100-x (Sb 2 Se 3 ) x chalcogenide glasses.

For a comparison, we employed IR VASE for a determination of refractive index of commercially available As 2 Se 3 glass (AMTIR-2, Amorphous Materials Inc.). The results obtained are given in Fig. 3 showing very good agreement between our data and data provided by producer. We point out also excellent consistency of As 2 Se 3 glass refractive index data (taking into account absolute precision of the methods) between IR VASE and prism coupling technique [22] showing maximal Δn of 0.004 at 3390 nm. Finally, Fig. 3 shows matching of refractive index determined by VASE and producer data for Ge 28 Sb 12 Se 60 bulk glass (IG5, Schott AG); the data again agree very well (Δn=0.005 at maximum in MIR). Although VASE does not allow reaching the accuracy obtained in the case of the prism method, the data analysis realized for two commercial selenide glasses validates the measurement method of refractive index by mean of infrared ellipsometry.

Conclusions

Bulk (GeSe 2 ) 100-x (Sb 2 Se 3 ) x chalcogenide glasses (x = 5-60) were synthesized in high purity.

Glass transition temperatures decrease with increasing antimony content, while densities show opposite tendency. This behavior is connected with the changes of glassy matrix connectivity. The fundamental absorption edge of the glasses shifts to lower energies with increasing content of antimony selenide. On the other hand, linear refractive index is strongly rising at the same moment (from 2.41 for (GeSe 2 ) 95 (Sb 2 Se 3 ) 5 glass to 3.01 for (GeSe 2 ) 40 (Sb 2 Se 3 ) 60 composition, at 1550 nm). We point out good accordance between refractive indices data extracted from spectroscopic ellipsometry and prism coupling technique. The main finding of this work is the applicability of infrared spectroscopic ellipsometry for precise determination of refractive index of chalcogenide glasses (as demonstrated for As 2 Se 3 and Ge 28 Sb 12 Se 60 glasses) in broad spectral range covering NIR and MIR regions.

Figure captions

Fig. 1 .

 1 Fig. 1. Spectral dependences of absorption coefficient of (GeSe 2 ) 100-x (Sb 2 Se 3 ) x glasses.

Fig. 2 .

 2 Fig.2. Refractive index dispersion of bulk (GeSe 2 ) 100-x (Sb 2 Se 3 ) x glasses obtained by ellipsometry covering UV-VIS-NIR-MIR spectral range: 300-2300 nm (full curves) and 1700-20000 nm (dashed curves). Note that full and dashed curves are overlapped in the region of 1700-2300 nm.

Fig. 3 .

 3 Fig. 3. Refractive index dispersion of bulk As 2 Se 3 (full curves) and Ge 28 Sb 12 Se 60 (dashed curves) glasses determined by ellipsometry covering UV-VIS-NIR-MIR spectral range (300-20000 nm). For a comparison, data published by producer of commercial As 2 Se 3 (AMTIR-2) and Ge 28 Sb 12 Se 60 (IG5) glass are given. Inset shows 900-14000 nm spectral range in more detail.

  ) 100-x (Sb 2 Se 3 ) x glasses, their cut-off wavelengths (at 10 cm -1 ), refractive indices determined by prism coupling technique (at 1311 and 1551 nm, ±1.10 -3 ) and VASE (1550 nm, ±1.10 -2 ), and optical band-gap values extracted from VASE data by Cody-Lorentz model (±0.01 eV).

  Figure 1 Click here to download high resolution image

Table captions Table 1 .

 captions1 Theoretical chemical composition of fabricated (GeSe 2

1 Table 1 .

 11

	Sample	Theoretical composition	Cut-off wavelength (nm)	n Prism coupling 1311 nm 1551 nm	n VASE 1550 nm	Eg (eV)
	x=5	Ge 30.6 Sb 3.2 Se 66.1	678	2.423	2.414	2.42	2.17
	x=10	Ge 28.1 Sb 6.3 Se 65.6	689	2.483	2.471	2.47	2.11
	x=20	Ge 23.5 Sb 11.8 Se 64.7	749	2.579	2.563	2.56	2.02
	x=30	Ge 19.4 Sb 16.7 Se 63.9	792	2.688	2.672	2.67	1.86
	x=40	Ge 15.8 Sb 21.1 Se 63.2	843	2.801	2.780	2.78	1.78
	x=50	Ge 12.5 Sb 25 Se 62.5	885	2.931	2.904	2.89	1.70
	x=60	Ge 9.5 Sb 28.6 Se 61.9	934	3.042	3.012	3.01	1.49

Table 1

 1

Acknowledgements

The Czech Science Foundation (Project No. 13-05082S) and Ministry of Education, Youth, and Sports of the Czech Republic (Project CZ.1.07/2.3.00/30.0058) and the CNRS PICS (Projet International de Cooperation Scientifique) program financially supported this work.


Hugo Chauvet

Olivier Devauchelle

François Métivier

Eric Lajeunesse

Angela Limare

Recirculation cells in a wide channel

come

Recirculation cells in a wide channel H. CHAUVET, 1, a) O. DEVAUCHELLE, 1 F. METIVIER, 1 E. LAJEUNESSE, 1 and A. LIMARE [START_REF] Forterre | Stability analysis of rapid granular chute flows: formation of longitudinal vortices[END_REF] Laboratoire de dynamique des fluides géologiques, Institut de physique du globe de Paris (IPGP), Paris, FRANCE.

(Dated:) Secondary flow cells are commonly observed in straight laboratory channels, where they are often associated with duct corners. Here, we present velocity measurements acquired with an Acoustic Doppler Current Profiler in a straight reach of the Seine river (France). We show that a remarkably regular series of stationary flow cells spans across the entire channel. They are arranged in pairs of counter-rotating vortices aligned with the primary flow. Their existence away from the river banks contradicts the usual interpretation of these secondary flow structures, which invokes the influence of boundaries. Based on these measurements, we use a depth-averaged model to evaluate the momentum transfer by these structures, and find that it is comparable with the classical turbulent transfer. They appear in rapid granular flows [START_REF] Forterre | Stability analysis of rapid granular chute flows: formation of longitudinal vortices[END_REF] , in Rayleigh-Bénard-Poiseuille flows 2 , in the superficial layers of sea and lakes [START_REF] Jeffreys | On the transverse circulation in streams[END_REF] , in straight tubes [START_REF] Znaien | Experimental and numerical investigations of flow structure and momentum transport in a turbulent buoyancy-driven flow inside a tilted tube[END_REF][START_REF] Hallez | Turbulence-induced secondary motion in a buoyancy-driven flow in a circular pipe[END_REF] , as well as in straight channels and natural rivers [START_REF] Rodriguez | Laboratory measurements of 3-d flow patterns and turbulence in straight open channel with rough bed[END_REF][START_REF] Blanckaert | Influence of shallowness, bank inclination and bank roughness on the variability of flow patterns and boundary shear stress due to secondary currents in straight open-channels[END_REF][START_REF] Albayrak | Secondary currents and corresponding surface velocity patterns in a turbulent open-channel flow over a rough bed[END_REF] .

In laboratory channels, a broad spectrum of instruments provide detailed measurements of the velocity field. They reveal secondary flow structures in the plane normal to the primary flow [START_REF] Rodriguez | Laboratory measurements of 3-d flow patterns and turbulence in straight open channel with rough bed[END_REF][START_REF] Blanckaert | Influence of shallowness, bank inclination and bank roughness on the variability of flow patterns and boundary shear stress due to secondary currents in straight open-channels[END_REF][START_REF] Shiono | Turbulent open-channel flows with variable depth across the channel[END_REF][START_REF] Nezu | Turbulence in Open-Channel Flows, IAHR Monographs[END_REF] , among which are the stationary helicoidal flows observed in straight channels (Prandtl refers to them as "secondary flow of the second kind " 11 ). They were originally associated with duct corners [START_REF] Nezu | Turbulence in Open-Channel Flows, IAHR Monographs[END_REF][START_REF] Brundrett | The production and diffusion of vorticity in duct flow[END_REF]13 , where they were attributed to the turbulence anisotropy [START_REF] Jeffreys | On the transverse circulation in streams[END_REF][START_REF] Brundrett | The production and diffusion of vorticity in duct flow[END_REF][14][START_REF] Gessner | The origin of secondary flow in turbulent flow along a corner[END_REF][START_REF] Demuren | Calculation of turbulence-driven secondary motion in ducts with arbitrary cross section[END_REF][START_REF] Husser | Direct numerical simulation of turbulent flow in a square duct[END_REF] . Recent experiments reported their existence in wide channels (aspect ratio of about 10), where they arrange themselves as pairs of counter-rotating vortices aligned with the stream direction [START_REF] Rodriguez | Laboratory measurements of 3-d flow patterns and turbulence in straight open channel with rough bed[END_REF][START_REF] Blanckaert | Influence of shallowness, bank inclination and bank roughness on the variability of flow patterns and boundary shear stress due to secondary currents in straight open-channels[END_REF][START_REF] Albayrak | Secondary currents and corresponding surface velocity patterns in a turbulent open-channel flow over a rough bed[END_REF] .

To our knowledge, secondary flow cells have never been measured in very wide channels -the only configuration which separates unambiguously recirculation cells from corner vortices. Alluvial rivers are typically wide and shallow, with an aspect ratio of a few tens [START_REF] Metivier | Alluvial landscape evolution: what do we know about metamorphosis of gravel bed meandering and braided streams[END_REF] . Therefore, a significant proportion of their flow should be free from the influence of the banks. Nonetheless, various indirect observations such as sediment streaks [START_REF] Nezu | Turbulence in Open-Channel Flows, IAHR Monographs[END_REF][START_REF] Ikeda | Self-formed straight channels in sandy beds[END_REF][START_REF] Colombini | Longitudinal streaks[END_REF] or fluctuations in the suspended particles concentration [START_REF] Vanoni | Experiments on the transportation of suspended sediments by water[END_REF] have been associated to secondary flows in alluvial streams. There is however no clear consensus on the physical origin of these patterns [START_REF] Nezu | Turbulence in Open-Channel Flows, IAHR Monographs[END_REF][START_REF] Ikeda | Self-formed straight channels in sandy beds[END_REF][START_REF] Colombini | Longitudinal streaks[END_REF] .

To identify directly secondary flow cells in rivers, we need accurate measurements of the velocity field. Unfortunately, because they require heavy infrastructures, laboratory measurement techniques are of delicate use in the field. Instead, many river studies involve Acoustic Doppler Current Profilers (ADCP). An ADCP emits a diverging bunch of ultrasonic beams from a single source. The three-dimensional velocity field is then reconstructed from the measured radial velocities, thus making it possible to record velocity profiles with a light measurement setup. Due to the divergence of the beams, the velocity reconstruction relies on the assumption that the flow is uniform in the horizontal plane. This configurational limitation generally precludes ADCPs from measuring flow structures smaller than the water depth [START_REF] Gargett | Observing turbulence with a modified acoustic doppler current profiler[END_REF] . Recirculation cells, however, are still accessible to ADCP measurements, since they are streamwise invariant and stationary. Here we use these properties to identify recirculation cells in a straight reach of the Seine river (France). At this location, the river is wide and deprived of bed sediment, thus limiting the possible causes of recirculation. Based on these measurements, we then evaluate the momentum transfer by secondary flow cells, and compare it to the classical turbulent transfer.

II. MEASUREMENTS A. Acoustic Doppler Current Profiler

Our ADCP (model RDI-Workhorse 1200 kHz) emits four divergent acoustic beams. Each beam measures the radial velocities u r every 25 cm through the water column. The beam angle with respect to the vertical direction is α = 20 • . Consequently, the distance 2l cor between opposite beams is proportional to the measurement depth h:

2l cor = 2h tan (α) ≃ 0.8h . (1)
Used conventionally, an ADCP computes the Cartesian components of the velocity field (u x , u y and u z ) from the four radial velocities (figure 1b). For instance, to get the streamwise velocity u x at a depth h, it uses the radial velocities from beam 3 and beam 4 (figure 1b).

The radial velocities u r,1 and u r,2 are decomposed as follows:

  u r,3 u r,4   =   u x,3 (-l cor ) u z,3 (-l cor ) -u x,4 (l cor ) u z,4 (l cor )   •   sin(α) cos(α)   , (2)
where u z,n and u x,n are respectively the vertical and the streamwise velocity in the n-th bins of each beam. This procedure relies on the assumption that the flow is uniform in the horizontal plane [START_REF] Gordon | Principles of operation a practical primer[END_REF] :

u i,3 (-l cor ) = u i,4 (l cor ). (3)
The streamwise velocity at deph h thus reads u x = (u r,3 -u r,4 )/2 sin(α). Similarly, beam 1 and beam 2 provide the transverse velocity u y 24 .

The uniformity assumption proves reasonable to measure integrated quantities, such as the river discharge, or large flow structures. However, due to the large tilt angle α of a beam, the correlation length l cor is of the order of the flow depth (equation 1, figure 1). Therefore, it cannot detect flow patterns smaller than the flow depth [START_REF] Gargett | Observing turbulence with a modified acoustic doppler current profiler[END_REF] . This precludes ADCPs from measuring the complete velocity field associated to secondary flow cells in a river.

In a straight channel however, we expect the average flow to be uniform in the streamwise direction, even in the presence of stationary helicoidal cells. The uniformity assumption therefore holds along the primary flow and, in principle, we can calculate the vertical velocity of secondary flow cells from the two acoustic beams aligned with the flow (beam 3 and beam 4 in figure 1a) [START_REF] Stacey | Measurements of reynolds stress profiles in unstratified tidal flow[END_REF][START_REF] Nystrom | Evaluation of mean velocity and turbulence measurements with ADCPs[END_REF] .

Secondary flow cells are weak stationary eddies in a highly turbulent flow. We thus need to average the velocity field over time to reveal them. In addition, since the uniformity assumption holds only on average, our experimental setup can only measure the average velocity field in the vertical streamwise plane (x, z). We define the average velocity as

u i = 1 T T 0 u i dt , (4)
where i denotes the coordinate x or z, and T is a time period long enough to average out the unsteady turbulent flow (T ≈ 10 min). Fortunately, u x and u z suffice to identify secondary flow cells, and to estimate their influence on the flow.

B. Field site

We collected our measurements midway of a 3 km-long straight reach of the Seine river (France). The river flows above bedrock between two vertical banks (width: 148 m, depth: about 6.2 ± 0.3 m depending on the discharge, figure 2a). The ADCP was mounted on a raft moored to a suspension footbridge which does not perturb the flow (figure 2b). The horizontal velocity of the raft was recorded by an echo-sounder, and the measurements are corrected for it. The raft position was measured with an auto-tracking theodolite with a precision of about a centimeter. We also used the ADCP goniometer to correct the measurements for the angular motion of the raft.

A. Streamwise velocity

To visualize the primary flow, we define the depth-averaged streamwise velocity as

U x (y) = 1 H H 0 u x dz , ( 5
)
where H is the water depth. The primary velocity appears roughly constant around the center of the channel, and slows down near the banks (figure 3a). As a result, the velocity profile is curved over about two thirds of the river width. Since the water depth is virtually constant across the river, this profile suggests that the flow transfers streamwise momentum to the banks.

- The vertical velocity profile is also affected by the banks. At the center of the channel, the velocity profile resemble the classical logarithmic law of the wall (figure 3b). However, as we approach the banks, the vertical profile deforms: the point of maximum velocity detaches from the water surface and migrates downwards, until it reaches approximately half the flow depth (figure 3c-d).

This observation is reminiscent of laboratory experiments showing a velocity maximum at about z = 0.6H near the banks [START_REF] Nezu | Turbulence in Open-Channel Flows, IAHR Monographs[END_REF][START_REF] Gibson | On the depression of the filament of maximum velocity in a stream flowing through an open channel[END_REF] . This departure from the logarithmic profile is usually interpreted as a signature of secondary flow cells near duct corners [START_REF] Blanckaert | Influence of shallowness, bank inclination and bank roughness on the variability of flow patterns and boundary shear stress due to secondary currents in straight open-channels[END_REF][START_REF] Nezu | Turbulence in Open-Channel Flows, IAHR Monographs[END_REF][START_REF] Prandtl | Essentials of fluid dynamics: With applications to hydraulics aeronautics, meteorology, and other subjects[END_REF] : these cells create an upward velocity along the side wall, thus bringing slow fluid up to the surface. The center of the cell, hardly impacted by the cell velocity, flows faster than the surrounding fluid. The same flow cell would have a much weaker effect far away from the banks, since the slow water must be transported from the bottom to the surface, thus explaining why the profiles near the center and near the bank differ.

Similarly, the weak spatial oscillation of the streamwise velocity across the river suggests that periodic upwelling brings slow water from the bottom to the surface (figure 3a).

B. Vertical velocity

The streamwise velocity profiles suggests that secondary flow cells exist across the entire channel. We now turn our attention to the vertical velocity, which provides direct insight on such flow structures (figure 4).

To reduce the scatter, the data set is divided into 2 m-wide bins. The vertical velocity is then averaged over each bin, and the standard deviation within a single bin provides an estimate of the error on the mean. The mean vertical velocity shows periodic oscillations across the river, with an amplitude significantly larger than the measurement error. As expected, the average vertical velocity across the channel vanishes, indicating that the vertical motion of the raft does not perturb significantly the measurements.

Despite considerable scatter, both the amplitude and the wavelength of this oscillation are roughly constant over the right-hand half of the river (the duration of the high-flow conditions did not allow measurements across the entire channel). They also show a reasonable consistency between data sets collected one year apart.

Fitting a cosine to the data leads to a wavelength of about 12 ± 0.5 m, and an amplitude of about 0.4 ± 0.1 cm s -1 . The Fourier spectrum of the same data set shows a maximum at a wavelength of about 12 ± 0.5 m (figure 5). Based on this observation, and in accordance with previous experiments 7, [START_REF] Albayrak | Secondary currents and corresponding surface velocity patterns in a turbulent open-channel flow over a rough bed[END_REF] , we suggest that the flow generates recirculation cells across the entire channel. The wavelength of the velocity oscillation corresponds to a series of counterrotating vortices with a diameter of about the flow depth (figure 4b and figure 6). The amplitude of the vertical oscillation is about 0.3% of the streamwise velocity, in accordance with previous findings in laboratory flumes [START_REF] Blanckaert | Influence of shallowness, bank inclination and bank roughness on the variability of flow patterns and boundary shear stress due to secondary currents in straight open-channels[END_REF] .

IV. MOMENTUM TRANSFER BY RECIRCULATION CELLS

In this section, our intention is to quantify the contribution of secondary flow cells in transporting momentum across the flow. To do so, we use a heuristic model to include the advection by secondary flow cells in the streamwise momentum balance.

A. Momentum balance

In steady state, the Reynolds-averaged Navier-Stokes equation for streamwise momentum reads

ρ ∂ u j u x ∂x j = ∂τ x,j ∂x j + ρgS, (6)
where g, S and ρ are the acceleration of gravity, the river slope and the density of water respectively. The repeated j index implies summation over the transverse y and vertical z coordinates. The brackets denote ensemble average and τ x,j is the stress tensor, including the Reynolds components.

In a rectangular channel with a steady free surface, the integration of equation ( 6) over depth leads to ∂ ∂y

H 0 (ρ u x u y -τ y,x ) dz = -τ bottom + ρgHS. (7)
In the above equation, the left-hand term represents the cross-stream momentum transfer by the conjugate action of the secondary flow and turbulence [START_REF] Shiono | Turbulent open-channel flows with variable depth across the channel[END_REF] . Even though the two integrands ρ u x u y and τ y,x can be isolated formally, they are physically interdependent.

The bottom shear stress τ bottom is often empirically approximated by a friction law:

τ bottom ≈ ρC f U 2 , ( 8
)
where U is the mean velocity and C f is a friction coefficient.

B. Turbulent transfer of momentum

We first consider the contribution of turbulence to the cross-stream transfer of momentum. Discarding temporarily the influence of the secondary flow, we use an empirical eddy viscosity ν t to get a rough estimate of the turbulent transfer:

τ x,y ≈ ρν t ∂ u x ∂y . (9)
Based on laboratory and field measurements, the vertical profile of eddy viscosity in an open channel is often approximated by a parabola 10,28-30 :

ν t ≡ κ u ⋆ z 1 - z H , ( 10
)
where κ is the von Kármán constant and u ⋆ is the friction velocity defined as ρu 2 ⋆ = τ bottom . In an infinitely wide channel (that is, without any flux of momentum across the stream), the velocity profile resulting from equation ( 6) with a turbulent viscosity closure ( 10) is logarithmic:

u x = u ⋆ κ ln z z 0 . (11)
where z 0 is the roughness length.

In a bounded channel, on the other hand, momentum is transferred from the flow bulk to the banks and, strictly speaking, equation (11) does not hold. However, if the channel aspect ratio is sufficiently large, we may assume that the vertical velocity profile remains logarithmic, in tune with the classical shallow-water approximation:

u x = U (y) C f κ ln z z 0 , (12)
where U (y) is the depth-averaged velocity.

Based on equations ( 9) and ( 12) , the cross-stream flux of momentum associated to the Reynolds stress tensor then reads

H 0 τ x,y dz ≈ B C f H 2 U ∂U ∂y (13)
where B depends on the ratio of the water depth to the roughness length only: Using the above expression to model the cross-stream momentum transfer, the momentum balance (7) reads

B = 1 6 ln H z 0 - 5 36 . (14) z y x
gHS = C f U 2 - BC f H 2 2 ∂ 2 U 2 ∂y 2 , (15)
To evaluate B from equation ( 14), we need an estimate of the roughness length z 0 . Fitting a logarithm on several vertical velocity profiles near the center of the Seine river (figure 3b), we find an average value of z 0 ≈ 2 ± 1 mm (this value is comparable with prior field measurements [START_REF] García | Sedimentation engineering: processes, measurements, modeling, and practice[END_REF] ), and therefore B ≈ 1.2 ± 0.1.

C. Influence of the the secondary flow on the momentum transfer

We now consider the influence of the secondary flow on the cross-stream transfer of momentum. Based on both the fluctuation of the vertical velocity across the stream (figure 4), and mass balance, we expect the recirculation cells to be arranged in a series of counterrotating eddies, with a diameter of about the water depth and a velocity U t (figure 6). Such a flow configuration, by itself, cannot transfer momentum beyond one cell diameter, as it is compartmented by vertical planes of vanishing cross-stream velocity. However, turbulent diffusion can take over the momentum transfer across these vertical planes, where the momentum gradient is enhanced by the recirculation cells. The mixing due to counter-rotating eddies could thus increase the diffusion of momentum across the stream [START_REF] Audoly | Réaction diffusion en écoulement stationnaire rapide[END_REF] .

Assuming that the momentum transfer is limited by convection within the cells, we expect it to scale like ρH U t (U --U + ) where U -and U + stand for the streamwise velocity left and right of a cell, respectively (figure 6). If the streamwise velocity varies on scales much larger than the water depth, we can approximate the velocity difference across a cell by a gradient, namely U --U + ≈ H∂U/∂y. This approximation is valid only in channels much larger than their depth, which is usually true for alluvial rivers [START_REF] Metivier | Alluvial landscape evolution: what do we know about metamorphosis of gravel bed meandering and braided streams[END_REF] .

Finally, based on prior experimental observations [START_REF] Blanckaert | Influence of shallowness, bank inclination and bank roughness on the variability of flow patterns and boundary shear stress due to secondary currents in straight open-channels[END_REF][START_REF] Nezu | Turbulence in Open-Channel Flows, IAHR Monographs[END_REF] and on the measurements presented in section III, we assume that the velocity of the secondary flow scales like that of the primary flow:

U t ≈ C d U (y) , (16)
where C d ≈ 0.003 ± 0.001. The momentum balance then reads

gHS = C f U 2 - C d H 2 2 ∂ 2 U 2 ∂y 2 . ( 17
)
This expression is similar to equation ( 15), except the dimensionless coefficient C d is substituted for BC f .

D. Comparison with field data

We now proceed to compare the estimates of the cross-stream flux of momentum to our field measurements in the Seine river. Equations ( 15) and (17), with the requirement that the velocity vanishes at the banks, share the following analytical solution:

U x = U ∞ 1 - cosh (y/L t ) cosh (w/(2L t )) , (18)
where

U ∞ = g H S C f . (19)
The diffusion length L t determines the inflection of the velocity profile near the banks. Its mathematical expression depends on the cross-stream diffusion model:

L t = H B 2
(no secondary flow, equation ( 15)) (20)

L t = H C d 2C f (with secondary flow, equation (17)). (21)
We evaluate the parameters U ∞ ≈ 1.2 m s -1 and L t ≈ 14 m by fitting equation (18) to the data (figure 7). The resulting shape of the velocity profile accords reasonably with the data. Since the the slope of the Seine river in Paris is about 10 -4 , we find C f ≈ 0.0042 ± 0.0001, which is consistent with classical empirical formulas [START_REF] Chezy | Formule pour trouver la vitesse de l'eau conduite dans une rigole donnée[END_REF][START_REF] Yalin | Fluvial Processes[END_REF] . The estimation of U ∞ , and therefore of C f , is virtually independent from the value of the diffusion length L t .

Based on our estimate of C f and on relations (20) and ( 21), we can calculate the diffusion length for the two diffusion models:

L t ≈ 4.8 ± 0.1 m (no secondary flow) (22)
L t ≈ 3.7 ± 0.6 m (with secondary flow).

The estimated diffusion length is similar for both models, and about three times smaller than the fitted value. Therefore, the two diffusion models not only share the same mathematical expression (equations ( 15) and ( 17)), but also involve similar constants. We thus expect the turbulent transfer of momentum to be comparable with the transfer induced by the secondary flow.

The similarity between the two transfer mechanisms implies (i) that recirculation cells cannot be neglected in the momentum balance, and (ii) that the momentum transfer by the recirculation cells cannot be distinguished from the turbulent transfer in the shallow-water framework presented here.

V. DISCUSSION AND CONCLUSION

Due to the diverging configuration of its acoustic beams, an ADCP generally cannot detect flow structures smaller than the flow depth. However, if the measurements are averaged over a sufficient period of time, this instrument can measure stationary and spatially extended flow structures.

We can unambiguously identify secondary flow cells in the vertical velocity field of the Seine River. These stationary vortices extend across the entire channel and rotate at about 0.3 % of the streamwise velocity, in accord with prior laboratory observations [START_REF] Rodriguez | Laboratory measurements of 3-d flow patterns and turbulence in straight open channel with rough bed[END_REF][START_REF] Blanckaert | Influence of shallowness, bank inclination and bank roughness on the variability of flow patterns and boundary shear stress due to secondary currents in straight open-channels[END_REF] .

An order-of-magnitude analysis suggests that the momentum transfer by these secondary In rivers, the momentum distribution determines stress on the bed, and therefore controls sediment transport and bed erosion. In addition, secondary cells are likely to enhance the horizontal mixing of passive markers, such as suspended particles, temperature or solutes [START_REF] Bouchez | Turbulent mixing in the Amazon River: The isotopic memory of confluences[END_REF] .

Finally, beyond open-channel flows, similar recirculation cells appear in a variety of turbulent sheared flows, such as longitudinal streaks in boundary layers [START_REF] Wanmin | Turbulent boundary-layer flow over fixed aerodynamically rough two-dimensional sinusoidal waves[END_REF][START_REF] Aihara | Formation of longitudinal vortices in the sublayer due to boundary-layer turbulence[END_REF] , Langmuir circulations below the wind-blown surface of lakes and oceans [START_REF] Langmuir | Surface motion of water induced by wind[END_REF] , or even granular flows 39 . To this day, we do not know how many of these phenomena share a common physical origin, if any.

(Received 21

 21 August 2013; accepted 19 December 2013; published online 23 January 2014) [http://dx.doi.org/10.1063/1.4862442] a) Electronic mail: chauvet@ipgp.fr; http://www.ipgp.fr/∼chauvet I. INTRODUCTION Turbulent shear flows often generate stationary vortices aligned with the primary flow.

FIG. 1 .

 1 FIG. 1. (a) Configuration of the ADCP beams, oriented with respect to the primary flow. (b) Configuration of beams 4 and 3, in the vertical streamwise plane. The tilt angle of the beam is α = 20 • and its opening angle is φ = 30 • . The radial velocity is measured at depth h. The associated correlation length is l cor .

FIG. 2 .

 2 FIG. 2. (a) Aerial picture of the field site (source: Institut Géographique National ). Measurements are made from the footbridge. (b) Set-up of the measurement raft.

  FIG. 3. Streamwise velocity measurements (m.s -1 ) in the Seine river during winter high flow. The discharge is 900 m 3 /s. (a) Depth-averaged velocity U x from moving measurements (blue) and from static measurements (yellow). (b), (c), (d) Time-averaged (≃ 10 minutes) depth profiles (u x ) with fitted logarithmic profiles (red curve).

FIG. 4 .

 4 FIG. 4. (a) Time-averaged vertical velocity u z in a transect, at half the flow depth. The shaded envelopes represent the standard error on the mean. The red curve is a cosine fitted to the data. (b) Topography of the river cross-section at the measurement location (scales are preserved). The distance y is measured from the middle of the channel.

FIG. 5 .

 5 FIG. 5. Energy spectrum of the vertical velocity field, calculated from the data set of figure 4.

FIG. 6 .

 6 FIG.6. Schematic representation of the recirculation cells.

FIG. 7 .

 7 FIG.7. Comparison of the depth-averaged momentum balance (red curve) with the ADCP measurements (yellow dots). Equation (18) is represented with U ∞ ≈ 1.2 m s -1 and L t ≈ 14 m fitted to the data.

  flow cells could compare with the intensity of the turbulent stress. If confirmed, this would indicate that the secondary flow influences significantly the primary flow in open channels. This is a strong incentive for further investigations. For instance, detailed flow measurements in laboratory experiments could reveal the streamwise extension of the flow cells, their stability in time and how they depend on the Reynolds number. Experiments would also guide us towards the mechanism which generates these structures.

B. Galletti and A. Bottaro, "Large-scale secondary structures in duct flow," Journal of

ACKNOWLEDGMENTS

The authors would like to thank the Navigation Service of the Seine river, represented by X. Fromageau, who granted access to the field site. A. Vierra built the measurement set-up.

We also thank M. Rossi, P.Y. Lagrée and F. Moisy for fruitful discussions.