Andrea Girón

Hugo Saiz

Flora S Bacelar

Roberto F S Andrade

Jesus Gómez-Gardeñes

Jesús Gómez-Gardeñes

Synchronization unveils the organization of ecological networks with positive and negative interactions

Keywords: numbers: 02.50.-r, 87.23.-n, 89.65.-s, 89.75.-k

come   L'archive ouverte pluridisciplinaire

In the last decade network science has provided the ideal benchmark to encode, analyze and understand the complex relationships that are stablished in large scale systems of nature as disparate as the Internet or the brain. As databases become more abundant and complete we face the challenge of analyzing networks containing multiple types of connections. Here we tackle the case of signed networks, where the interactions can be either positive or negative. On one hand we illustrate how synchronization processes capture the organization of this kind of graphs into a set of modules interacting negatively among them. With this in mind, we apply the synchronization benchmark to real systems. In particular, we characterize the ecological balance between facilitation and competition in plant communities and analyze the role of species in their organization.

I. INTRODUCTION

Synchronization is perhaps the most paradigmatic example of collective behavior as it is recurrently found at different levels of complexity [START_REF] Pikovsky | Synchronization: a Universal Concept in Nonlinear Sciences[END_REF][START_REF] Manrubia | Emergence of dynamical order: synchronization phenomena in complex systems[END_REF][START_REF] Boccaletti | The Synchronized Dynamics of Complex Systems[END_REF]. In fact, the emergence of spontaneous synchronization in systems of coupled dynamical units is at the core of many coordinated tasks, from cognitive processes in the brain, to the unfolding of collective behaviors in social systems [START_REF] Watts | Collective dynamics of smallworld networks[END_REF][START_REF] Strogatz | Exploring complex networks[END_REF][START_REF] Bullmore | The economy of brain network organization[END_REF]. In the last decade network theory has unveiled that the topology of the interactions in a complex system has important effects on the development of collective behaviors [START_REF] Strogatz | Exploring complex networks[END_REF][START_REF] Boccaletti | Complex Networks: Structure and Dynamics[END_REF]. Following this direction, the study of synchronization in networks has attracted a lot of attention [START_REF] Arenas | Synchronization in complex networks[END_REF] to shed light on the role that the network structure plays on the emergence of synchronized states [START_REF] Nishikawa | Heterogeneity in oscillator networks: are smaller worlds easier to synchronize?[END_REF][START_REF] Moreno | Synchronization of Kuramoto oscillators in scale-free networks[END_REF][START_REF] Chavez | Synchronization is enhanced in weighted complex networks[END_REF][START_REF] Arenas | Synchronization reveals topological scales in complex networks[END_REF][START_REF] Zhou | Universality in the Synchronization of Weighted Random Networks[END_REF][START_REF] Gómez-Gardeñes | Paths to synchronization in complex networks[END_REF][START_REF] Lodato | Synchronization Properties of Network Motifs[END_REF][START_REF] Gómez-Gardeñes | Explosive Synchronization Transitions in Scale-Free Networks[END_REF][START_REF] Gambuzza | Analysis of remote synchronization in complex networks[END_REF].

The typical setting of the former works consists on associating a dynamical system to each node, whereas the couplings between pairs of dynamical units are mediated by the links of the network. However, other studies have also covered adaptive networks [START_REF] Gross | Adaptive coevolutionary networks: a review[END_REF] whose structure is shaped by the microscopic synchronization patterns [START_REF] Zhou | Dynamical Weights and Enhanced Synchronization in Adaptive Complex Networks[END_REF][START_REF] Sorrentino | Adaptive Synchronization of Dynamics on Evolving Complex Networks[END_REF][START_REF] Aoki | Co-evolution of Phases and Connection Strengths in a Network of Phase Oscillators[END_REF][START_REF] Gutierrez | Emerging meso-and macroscales from synchronization of adaptive networks[END_REF] or systems of mobile oscillators moving in a continuous space [START_REF] Buscarino | Dynamical network interactions in distributed control of robots[END_REF][START_REF] Frasca | Synchronization of Moving Chaotic Agents[END_REF][START_REF] Fujiwara | Synchronization in networks of mobile oscillators[END_REF][START_REF] Prignano | Tuning Synchronization of Integrate-and-Fire Oscillators through Mobility[END_REF][START_REF] Gómez-Gardeñes | Motion-induced synchronization in metapopulations of mobile agents[END_REF]. With relatively few exceptions [START_REF] Daido | Quasientrainment and slow relaxation in a population of oscillators with random and frustrated interactions[END_REF][START_REF] Stiller | Dynamics of nonlinear oscillators with random interactions[END_REF][START_REF] Daido | Algebraic relaxation of an order parameter in randomly coupled limit-cycle oscillators[END_REF][START_REF] Leyva | Sparse repulsive coupling enhances synchronization in complex networks[END_REF][START_REF] Nishikawa | Network synchronization landscape reveals compensatory structures, quantization, and the positive effect of negative interactions[END_REF][START_REF] Hong | Model of Coupled Oscillators with Positive and Negative Coupling Parameters: An Example of Conformist and Contrarian Oscillators[END_REF][START_REF] Hong | Conformists and contrarians in a Kuramoto model with identical natural frequencies[END_REF][START_REF] Hong | Mean-field behavior in coupled oscillators with attractive and repulsive interactions[END_REF][START_REF] Anderson | Multiscale dynamics in communities of phase oscillators[END_REF], the hypothesis of these works is that the interactions between units are positive, so that the existence of a link between two coupled dynamical systems implies that they are prone to synchronize as the interaction between them is increased.

In this work, we address the scenario in which positive and negative interactions between nodes coexist in the same network. In this way, the increase of the interaction strength causes the attraction among those dynamical units interacting via positive links and a repulsive effect between those connected through negative edges. As shown recently by Anderson and co-workers [START_REF] Anderson | Multiscale dynamics in communities of phase oscillators[END_REF], the stable equilibrium reached consists in a dynamical partition of the network in which nodes sharing positive connections synchronize together while they avoid being dynamically close to those with whom negative interactions are at work. Here, our goal is to use this partition to extract information about the organization of networks for which positive and negative relationships between nodes interplay. In particular we will analyze plant communities where usually both facilitation (positive interaction) and competition (negative interaction) coexist [START_REF] Callaway | Competition and facilitation: a synthetic approach to interactions in plant communities[END_REF], playing a key role in the structure and organization of these communities [START_REF] Aguiar | Patch structure, dynamics and implications for the functioning of arid ecosystems[END_REF][START_REF] Bacelar | Exploring the tug of war between positive and negative interac-tions among savanna trees: Competition, dispersal, and protection from fire[END_REF][START_REF] Saiz | Changes in semi-arid plant species associations along a livestock grazing gradient[END_REF][START_REF] Tirado | Community structure and positive interactions in constraining environments[END_REF]. We will show that synchronization provides a efficient way to unveil how species sharing facilitative interactions group into blocks, allowing the identification of the species that play a major role in the structure of the plant community.

A c c e p t e d m a n u s c r i p t

The article is organized as follows. First in Sec. II and Sec. III we review the main results about the Kuramoto model of coupled phase oscillators and its studies in signed networks respectively. Then, in Sec. IV we analyze a toy model in which two networked systems of attractively coupled oscillators interact repulsively. This allow us to tackle the characterization of real plant communities in Sec. V exploring the emergence of facilitation groups and the competitive interactions established among them. Finally, in Sec. VI we round off the article by drawing the concluding remarks.

II. THE KURAMOTO MODEL

In the following, we rely on the paradigmatic Kuramoto model [START_REF] Kuramoto | Self-entrainment of a population of coupled nonlinear oscillators[END_REF][START_REF] Strogatz | to Crawford: exploring the onset of synchronization in populations of coupled oscillators[END_REF][START_REF] Acebrón | The Kuramoto model: A simple paradigm for synchronization phenomena[END_REF] (KM) of coupled limit circle oscillators as it is the most used dynamical framework for the study of synchronization phenomena in complex networks [START_REF] Arenas | Synchronization in complex networks[END_REF]. In a complex network of size N , the Kuramoto model considers that each node i is characterized by a phase θ i ∈ (0, 2π] whose time evolution is given by:

θi = ω i + λ N ∑ j=1 A ij sin(θ j -θ i ) .
(

) 1
The above equation thus describes a set of phase-oscillators coupled with uniform strength λ as dictated by the adjacency matrix of the network A, whose terms are A ij = A ji = 1 if nodes i and j are connected and A ij = 0 otherwise. Note that, in the original Kuramoto model [START_REF] Kuramoto | Self-entrainment of a population of coupled nonlinear oscillators[END_REF], each node i interacts with all the others nodes, so that

A ij = (1 -δ ij ).
The natural frequencies of the oscillators, {ω i }, are, in principle, different and they are assigned following a frequency distribution g(ω), that is usually considered to be uni-modal and even around the mean frequency Ω of the population.

In the uncoupled regime (λ = 0) each node i describes limit-cycle oscillations with characteristic frequency ω i . However, by increasing the coupling strength λ, the network undergoes a phase transition at some critical value λ c that, for an arbitrary network topology described by an adjacency matrix A, is given by λ c = 2/[πg(Ω)Λ(A)] [START_REF] Restrepo | Onset of synchronization in large networks of coupled oscillators[END_REF] where g(Ω) is the density of oscillators with natural frequency equal to the average one, Ω, and Λ(A) is the maximum eigenvalue of the adjacency matrix.

To monitor the synchronization transition Kuramoto introduced the following order parameter:

r(t)e iΨ(t) = 1 N N ∑ j=1 e iθj(t) , (2)
which is the averaged sum of the unit vectors associated to the phases of each oscillator in the complex plane. The modulus of the resulting complex number is the order parameter, r(t) ∈ [0, 1], that measures the coherence of the collective motion of the oscillators. In this way, r = 1 when the network is fully synchronized whereas r = 0 stands for the incoherent solution. In addition, the value of Ψ(t) accounts for the average phase of the ensemble of oscillators.

III. THE KURAMOTO MODEL IN SIGNED NETWORKS

The extension of the KM, Eq. ( 1), to address the possibility of positive and negative interactions was initially proposed by Daido [28] more than twenty years ago. There he considered an all-to-all Kuramoto model with random interaction weights (or equivalently random entries A ij ) including positive and negative ones. From this study a debate about the glassy behavior of the system was opened and its solution still remains unclear [START_REF] Stiller | Dynamics of nonlinear oscillators with random interactions[END_REF][START_REF] Daido | Algebraic relaxation of an order parameter in randomly coupled limit-cycle oscillators[END_REF]. Recently, the combination of positive and negative interactions in a network of Kuramoto oscillators has recovered the attention of researchers [START_REF] Nishikawa | Network synchronization landscape reveals compensatory structures, quantization, and the positive effect of negative interactions[END_REF][START_REF] Hong | Model of Coupled Oscillators with Positive and Negative Coupling Parameters: An Example of Conformist and Contrarian Oscillators[END_REF][START_REF] Hong | Conformists and contrarians in a Kuramoto model with identical natural frequencies[END_REF][START_REF] Hong | Mean-field behavior in coupled oscillators with attractive and repulsive interactions[END_REF]. For instance, in [START_REF] Nishikawa | Network synchronization landscape reveals compensatory structures, quantization, and the positive effect of negative interactions[END_REF] the authors show that the addition of negative links to networks can enhance the stability of the synchronized state, while Hong and Strogatz [START_REF] Hong | Model of Coupled Oscillators with Positive and Negative Coupling Parameters: An Example of Conformist and Contrarian Oscillators[END_REF][START_REF] Hong | Conformists and contrarians in a Kuramoto model with identical natural frequencies[END_REF][START_REF] Hong | Mean-field behavior in coupled oscillators with attractive and repulsive interactions[END_REF] studied two different coarse-grained versions of the Daido model to shed light on the aforementioned debate.

Here we will revisit the KM in the following way. Instead of working with weighted adjacency matrices we allow the weights of the links to take one out of three values. Namely, A ij = A ji = 1 if the interaction between i and j is positive, A ij = A ji = -1 if the link between i and j is negative and A ij = A ji = 0 if no interaction between i and j is present. This is equivalent to work with two adjacency matrices, one encoding all the positive interactions, A + , and another one, A -, accounting for the negative ones. In this way we have A = A + -A -, and the KM, Eq. ( 1), can be rewritten as:

θi = ω i + λ N ∑ j=1 [ A + ij sin(θ j -θ i ) + A - ij sin(θ i -θ j ) ] . (3)
The above equation is the most simple generalization of the KM in a complex network to the case of signed graphs. However, it can be applied to many real situations in which positive and negative interactions coexist within the same system. A familiar case is that of social networks, in which the acquaintances of an individual can be classified as friends and foes. In this latter case, the focus has been the test of social balance theory [START_REF] Heider | Attitudes and cognitive organization[END_REF][START_REF] Cartwright | Structure balance: A generalization of Heider's theory[END_REF] in large scale social networks [START_REF] Szell | Multirelational organization of large-scale social networks in an online world[END_REF] where complex structures mixing positive and negative ties appear.

As anticipated, here we will focus on a different real world example where the interplay between negative and positive interactions appears naturally: ecological networks. We rely on recent data sets about the interaction between vegetal species [START_REF] Saiz | Effect of livestock grazing in the partitions of a semiarid plant-plant spatial signed network[END_REF][START_REF] Saiz | Plant-plant spatial association networks in gypsophilous communities: the influence of aridity and grazing and the role of gypsophytes in its structure[END_REF] in which the nature of the interactions between the different species was derived by analyzing their local spatial association patterns. Although an approximation to real biotic interactions, spatial association between plant species is considered as a valid proxy for studies including all species in a community [START_REF] Bacelar | Exploring the tug of war between positive and negative interac-tions among savanna trees: Competition, dispersal, and protection from fire[END_REF][START_REF] Cavieres | Nurse effect of Bolax gummifera cushion plants in the alpine vegetation of the Chilean Patagonian Andes[END_REF][START_REF] Cavieres | Positive interactions between alpine plant species and the nurse cushion plant Laretia acaulis do not increase with elevation in the Andes of central Chile[END_REF]. These interactions can be of three types: neutral (no interaction, species associate at random in space), facilitation (positive interaction, species cooccur in space more than expected) or competition (negative interaction, species segregate in space more than expected). In this way, the two matrices A + and A -can be constructed straightforward.

A c c e p t e d m a n u s c r i p t IV. A SIMPLE SYSTEM: TWO INTERCONNECTED NETWORKS

The KM has been used to unveil the modular organization of unsigned networks [START_REF] Arenas | Synchronization reveals topological scales in complex networks[END_REF][START_REF] Gómez-Gardeñes | Synchronizability determined by coupling strengths and topology on Complex Networks[END_REF][START_REF] Boccaletti | Detecting complex network modularity by dynamical clustering[END_REF][START_REF] Almendral | Dynamics of overlapping structures in modular networks[END_REF][START_REF] Sendiña-Nadal | Unveiling Protein Functions through the Dynamics of the Interaction Network[END_REF] relying on the observation that synchronization is first attained among those nodes belonging to the same module, or community, and then spreads to the rest of the network. This provides a fast an efficient way to detect the modular structure of the network. For signed networks, Anderson et al. [START_REF] Anderson | Multiscale dynamics in communities of phase oscillators[END_REF] analyzed the behavior of a system composed of M large groups of coupled Kumamoto oscillators, so that the interactions within members of the same group are positive whereas those between oscillators belonging to different groups are negative. To illustrate their results we consider here the simplest case of M = 2 (analyzed in [START_REF] Montbrió | Synchronization of two interacting populations of oscillators[END_REF][START_REF] Rad | Topological measure locating the effective crossover between segregation and integration in a modular network[END_REF][START_REF] Li | Synchronization Interfaces and Overlapping Communities in Complex Networks[END_REF] for unsigned networks) groups, see Fig. 1. Thus, as in [START_REF] Anderson | Multiscale dynamics in communities of phase oscillators[END_REF], we assume that there is a perfect partition of the system into two sub-populations that we know a priori. Our aim is to observe how synchronization patterns are organized so to gain insight for the application to real systems in which the sub-population structure, if any, is unknown.

As shown in Fig. 1 we consider a toy multilayer network [START_REF] Boccaletti | The structure and dynamics of multilayer networks[END_REF][START_REF] Kivela | Multilayer networks[END_REF] composed two fully-connected graphs of identical size N , while the probability that two nodes belonging to different networks are connected is equal to p. In this way, each network contains N (N -1)/2 positive links while, on average, the two networks are interconnected via pN 2 links. Thus, the system contains N 2 positive and pN 2 negative links. Considering the probability of connection p between two nodes of different populations, the KM for the system of two interconnected networks can be approximated as:

θα i = ω α i + λ N ∑ j=1 sin(θ α j -θ α i ) -pλ N ∑ j=1 sin(θ β j -θ α i ) , ( 4
)
where θ α i and ω α i are the phase of oscillator i (i = 1, ..., N ) in population α and its natural frequency respectively. Note that in Eq. ( 4) α = 1, 2 and β = 2, 1 respectively. In addition to the Kuramoto order parameter r it is useful to define a synchronization parameter associated to each population α. In analogy to Eq. ( 2) this order parameter reads:

r α (t)e iΨ α (t) = 1 N N ∑ j=1 e iθ α j (t) , (5)
where r α ∈ [0, 1] accounts for the degree of synchronization in population α whereas Ψ α defines the average phase of oscillators belonging to network α. With the help of Eq. ( 5) we can re-write Eq. ( 4) as:

θα i = ω α i + Kr α sin(Ψ α -θ α i ) -pKr β sin(Ψ β -θ α i ) , ( 6
)
where K = N λ, recovering the usual mean-field form of the KM [START_REF] Kuramoto | Self-entrainment of a population of coupled nonlinear oscillators[END_REF][START_REF] Strogatz | to Crawford: exploring the onset of synchronization in populations of coupled oscillators[END_REF][START_REF] Acebrón | The Kuramoto model: A simple paradigm for synchronization phenomena[END_REF]. Assuming, without loss of generality, a reference frame rotating with the mean frequency Ω of the system (or analogously considering Ω = 0) and whose position is set so that the mean phase of the oscillators in one of the networks is Ψ α = 0, the equation for the locked, θα i = 0, oscillators reads:

ω α i = Kr α sin(θ α i ) + pKr β sin(Ψ β -θ α i ) . ( 7
)
The above equation is only fulfilled for a set of oscillators in each of the populations having a natural frequency close enough to the average one Ω = 0. Due to the high symmetry of the two populations and under the assumption that the frequency distribution is the same for both networks, g α (ω) = g β (ω), we can can establish a map between a locked oscillator i in population α, with its counterpart, say i, in population β so that ω α i = ω β i . In this way, for each pair of locked oscillators we have:

sin(θ α i ) + p sin(θ β i ) = sin(θ β i -Ψ β ) + p sin(θ α i -Ψ β ) (8)
where we have considered, again invoking the symmetry of the system and that g α (ω) = g β (ω), that the degree of synchronization of the two populations is identical, r α = r β , for a given value of K. The self-consistent set of equations ( 8) has two solutions:

Ψ β = 0 (θ α i = θ β i ) and Ψ β = π (θ α i = θ β i -π).
It can be easily shown (see Appendix A for a simple argument and the Supplemental Material [62] for a proof) that the first solution (Ψ α = Ψ β = 0) is unstable whereas the second one (Ψ α = 0,Ψ β = π) is stable, pointing out that although synchronization can be attained within each of the populations they tend to be in anti-phase so that the degree of synchronization of the whole system r falls to zero.

In Fig. 2 we show this phenomenon by considering two networks of N = 500 nodes and an interconnection probability of p = 0.4. The natural frequencies are randomly assigned from an homogeneous (flat) distribution within the range

ω α i ∈ [-1/2, 1/2].
In the top panels, (a)-(c), we show the distribution of the phases in the unit circle of the complex plane for K = 0, K = 0.5 and K = 1. It is shown how synchronization within the groups, pointed by the concentration of the circles, is attained as K is increased. From these panels it is also observed that the synchrony groups are in anti-phase. Panel (d), in the bottom, shows the evolution of the degree of synchronization within the groups, r 1 and r 2 , and that of the whole system, r, as a function of K. It becomes clear that the anti-phase alignment of the synchrony groups observed above makes r decrease to 0 while the values r α increase for K > K c ≃ 0.45. Now we focus on the onset of synchronization K c for the two populations. By considering the stable solution of Eq. ( 8) in Eq. ( 6) we can express the evolution for the phase of an oscillator i in population α as:

θα i = ω α i -(1 + p)Kr α sin(θ α i ) . ( 9
)
The above equation is exactly the same to the original KM with the exception of the (1 + p) correction in the coupling strength. Thus, within a population α synchronization shows up at some critical coupling λ c = K c /N where

K c = K KM c (1 + p) = 2 πg(0)(1 + p) . ( 10
)
In Fig. 3 we show the degree of synchronization within a network, r α , as a function of the coupling constant K normalized to the corresponding value in the original KM, K KM c . Obviously, for p = 0 we have two independent fully connected networks, i.e., two systems obeying the original KM, and thus the onset is reached exactly at K KM c . As predicted from Eq. (10) as p increases, so that negative links between networks appear, the onset of synchrony is anticipated, reaching half of the original value, K KM c , when p = 1. Apart from the knowledge gained from the analysis of the toy multilayer network, the anticipation of the onset of synchronization within each population as p increases constitutes a apparently counterintuitive result: the more repulsive links . The agreement with the theoretical prediction for the synchronization onset, Eq. ( 10), is excellent. The interconnected system is the same as in Fig. 2.

we add the easier synchronization is attained between nodes sharing attractive interactions. Translated into a social context this implies that two friends tend to become more similar as the number of common foes increases. This result constitutes a nice theoretical example of the importance that negative interactions have on the collective behavior, here synchronization, of complex systems.

V. UNCOVERING THE STRUCTURE OF REAL SIGNED NETWORKS

Now we address the characterization of two real ecological networks in which negative and positive interactions interplay. These networks capture the relationships established between vegetal species by measuring the spatial associations among them [START_REF] Saiz | Effect of livestock grazing in the partitions of a semiarid plant-plant spatial signed network[END_REF][START_REF] Saiz | Plant-plant spatial association networks in gypsophilous communities: the influence of aridity and grazing and the role of gypsophytes in its structure[END_REF]. Spatial associations are derived by analyzing the deviations with respect to the expected co-ocurrence frequency when plants are randomly distributed in space. In this way, those pairs of species co-occurring more than expected are considered to interact positively (facilitation), whereas those pairs co-occurring less than expected are considered to interact negatively (competition). The two studied networks correspond to two typical Iberian semi-arid locations, Cabo de Gata (Almería) and Monegros (Aragón), being both placed in Spain. The main properties of both networks are collected in Table I. In both networks positive interactions dominate over negative ones, a typical property of plant communities under stressful environmental conditions where positive interactions are hypothesized to be more important than negative ones [START_REF] Brooker | Facilitation in plant communities: the past, the present, and the future[END_REF][START_REF] Soliveres | Moving forward on facilitation research: response to changing environments and effects on the diversity, functioning and evolution of plant com-munities[END_REF]. The main difference between these two networks lies on their different density of links, being the Monegros networks the most wired one. A possible explanation is the presence of high livestock grazing intensity in the plant community of Cabo de Gata (0.65 ind/ha, for an area where grazing carrying capacity is 0.39-0.57 ind/ha, [START_REF] Robles | Native forage shrub species in south-eastern Spain: forage species, forage phytomass, nutritive value and carrying capacity[END_REF]). Overgrazing is known to have a deep impact in the structure of semi-arid plant communities, particularly breaking the spatial organization of vegetation patches and randomizing the associations between plant species [START_REF] Milchunas | Quantitative effects of grazing on vegetation and soils over a global range of environments[END_REF][START_REF] Alados | Association between competition and facilitation processes and vegetation spatial patterns in alpha steppes[END_REF]. Thus, the network is expected to present less linkage density. Nevertheless, the differences between the two networks allow us to validate the partitioning method in real communities with different properties.

To analyze how competitive and facilitation relationships among species interplay and shape the interaction network we proceed as follows. We set all the initial phases of the oscillators as θ i (0) = 0 ∀i and assign nearly identical natural frequencies ω i = 0 + ξ (with ξ being a random variable homogeneously distributed in the interval [-10 -3 , 10 -3 ]). Then we start to increase adiabatically the value of the interaction coupling K by small increments δK and monitor how the dynamical equilibria reached for each value of K evolve from small to large values of the coupling strength. Note that, at variance with the system analyzed in Section 4, here we do not know the partition, if any, of the network so that our goal is to detect how groups show up from the synchronization clusters.

A. Cabo de Gata Ecosystem

In Fig. 4 we analyze the synchronization patterns of the first network (Cabo de Gata). The first insight is provided by the evolution of the degree of synchronization, r, of the network as a function of the coupling K (see Fig. 4.a). Given the initial conditions described above for K ≪ 1, the oscillators remain fully synchronized (r = 1) as a product of the very small deviation in their natural frequencies. However, as K increases there is a sudden drop of r due to the repulsive effect of negative interactions and finally, for K > 0.2, the equilibrium state remains unaltered displaying a degree of synchrony of r ≃ 0.89.

Once each dynamical equilibrium is reached for a value of K we look at the distribution of the steady phases {θ ⋆ i } of the oscillators in the unit circle. In particular, we measure the relative phase of an oscillator i with respect to that of oscillator 1 (the most connected node in the network):

∆θ i = θ ⋆ i -θ ⋆ 1 .
In Fig. 4.b we observe that as soon as the coupling strength K is large enough a stable organization of the network into clusters of oscillators having the same dynamical state shows up. According to Fig. 4.b we can distinguish 6 clusters, denoted from I to VI in the right vertical axis Fig. 4.b. At first sight it becomes clear that all the species in groups II to V are strongly repelled by cluster I since these latter clusters are placed around ∆θ i = -π. In fact, by looking at the composition of cluster I we notice that this is composed by only node 1 (the hub of the system) whose interactions are all repulsive. Thus, node 1 appears as a very competitive species which play a central role structuring the plant community.

Taking advantage of the partition of the network provided by the Kuramoto dynamics we can organize the nodes of the network as in Fig. 4.c so to understand the composition and role of the remaining clusters II to VI. First we notice that all the members of cluster V are either direct competitors of node 1 (cluster I) or neighbors of these competitors, suggesting a large group of species that associate among them in the spatial gaps left by species 1. This result agrees well with the distribution shown in Fig. 4.b since cluster V is placed just at ∆θ i = -π whereas clusters II, III and IV, although also containing competitive interactions with node 1, are displaced in such a a way that ∆θ i > -π for the nodes they contain. The reason for such displacement comes from the competition of cluster II with cluster VI (composed of nodes 2 and 3) which, in its turn, also competes with cluster I (node 1). In fact, the competition between clusters I, II and VI manifests in their respective disposition observed in Fig. 4.b where they appear in such way that nearly maximize their relative distance. The organization of these clusters reveals the existence of two other competitors (cluster VI) not as strong as node 1 but with an specific negative effect on the species of cluster II. In principle, if these clusters I, II and VI had no interactions with clusters IV and III they would appear placed at ∆θ i = 0, -2π/3 and -4π/3 respectively. However, cluster II has positive interactions with cluster III which, in its turn, interacts negatively with cluster I. Thus, as a result of the attraction between clusters II and III, the positions of the nodes in these clusters appear displaced from ∆θ i = -2π/3 and -π respectively. Finally, the displacement of nodes in cluster III from ∆θ i = -π originates a (smaller) displacement of nodes in cluster IV from ∆θ i = -π due to the facilitation links [START_REF] Aoki | Co-evolution of Phases and Connection Strengths in a Network of Phase Oscillators[END_REF][START_REF] Zhou | Universality in the Synchronization of Weighted Random Networks[END_REF] and [START_REF] Buscarino | Dynamical network interactions in distributed control of robots[END_REF][START_REF] Chavez | Synchronization is enhanced in weighted complex networks[END_REF].

B. Monegros Ecosystem

Now we proceed to analyze the second network corresponding to the Monegros location. As shown in Table I and discussed above, this system is richer than that of Cabo de Gata in terms of the number of interactions and, at variance with the first network, facilitation interactions appear significantly more frequently than competitive ones. The analysis is done in the same fashion as for the first network (see Fig. 5). In Fig. 5.a we show the evolution of the degree of synchronization of the system as a function of K. This provides us a first insight about the minimum coupling strength needed to reach an equilibrium in which the competition between negative and positive interactions is manifested. This happens for K > 0.25. The evolution with K of the phase difference of each node i with respect to the most connected node (again node 1), ∆θ i , is shown in Fig. 5.b and corroborates that for K < 0.25 there is a self-organization of the phases into groups that become stable for K > 0.25.

Fig. 5.b clearly shows the organization of the species (nodes) into 10 groups (named as I, II,..., X) that are shown in the representation of the network in Fig. 5.c. The first group is composed by node 1 and all the species that interact (directly or indirectly) with it via positive links. Thus, node 1 plays a key role structuring the community by isolating species in group I from the rest of system. In this case, node 1 is clearly a nurse species responsible of patch formation. These patches act as shelters against the harsh environmental conditions, increasing the chances of survival of the species in cluster I.

The remaining nodes in the network compete (directly or indirectly) with node 1, and represent the species that are established outside the patches including species 1. However, this second part of the system appears fragmented into 9 The nodes of these networks represent groups of species whereas a positive (negative) link between a pair of groups exists if, at least, one positive (negative) interaction between two species belonging to each group is present in the original network. Solid (dashed) links reprsent positive (negative) interactions.

groups due to the internal competition links between the members of these groups. As competitors of node 1, the members of group V appear, as shown in Fig. 5.b, placed at ∆θ i = -π, whereas groups II to X are displaced from this antagonist position. The reason for the displacement lies on the competitive interactions launched from groups X and II which are competitors of group I as well. For instance, group X is competitor of groups I, III, IV, V so that it automatically implies that, considering the position of groups I (∆θ i = 0) and V (∆θ i = -π), the positions of groups III and IV lie in the range ∆θ i ∈ (0, -π) as they are also competitors of group I and are positively linked to group V.

The role of groups III and V would be identical (as given from their respective interactions with groups I, V and X) if it were not for the competitive interactions that group III shares with groups VII and VIII that, in their turn, are positively connected to group V. These negative links become evident from the position of group III in Fig. 5.b (∆θ i > -π) whereas group IV, not having these competitive interactions, lies close to group V (∆θ i ≃ -π). The negative interactions of group III with groups VII and VIII automatically set their position at ∆θ i < -π while being close to group V with whom share facilitation links. Moreover, the positive interactions of groups VII and VIII with groups IX and VI respectively set the latter close to the former, being group VI the closest to ∆θ i = -π of the four due to its positive link to group V. On the other hand, group IX is the farthest one due to the negative interaction with group II (another competitor of group I) which, as a consequence, appears in the region ∆θ i > -π.

C. Unveiling the Organization of modules

With the partition obtained above one can easily analyze the organization that facilitation and competitive interactions provides to the system. In Fig. 6 we show the coarse-grained representation of the Cabo de Gata (a) and Monegros (b) networks where the nodes account for the detected groups and the interactions between them are positive or negative depending on the nature of the interactions between the nodes each group contains. This picture provides a valuable tool to identify groups of species which present the same interaction behavior in the system. For example, a group of species which share positive interactions represents a particular type of vegetation patch in the community (as group I in the Monegros ecosystem). Let us remark that none of the pairs of connected groups share mixed interactions, i.e., the links between groups are either positive or negative. Thus, the partitions obtained from synchronization clusters allow us to interpret unambiguously the relationships between the different groups of species in an ecosystem.

For the first network (see Fig. 6.a), it is easy to notice the critical role that group VI has on the organization of the network. In fact, by deleting group VI the role played by II, III and IV would be identical so that they would merge into one single group and, analogously, this new resulting group would have an analogous role to that of group V implying a new merging of groups that would turn the system into a network of two competing populations. This conclusion can be easily reached by looking directly to the corresponding network of interactions between species (see Fig. 4.c). However, the analysis of the Monegros network (see Fig. 5.c) becomes more difficult.

In ecological terms, the pattern of interactions between the groups in the Monegros ecosystem, shown in Fig. 6.b, suggests the existence of a complex competition pattern among the modules competing with module I, probably related to the particular life histories of the species composing them. Specifically, by considering groups II to X we observe that groups II and X are monospecific so that the two species (22 and 20 respectively) establish themselves outside the vegetation patches composed of species belonging to groups I and III to IX. However, from Fig. 6.b it becomes clear that vegetation patches containing species from groups III to IX can be composed of species belonging to group V together with some species belonging to groups III, IV, and VI-IX. Specifically, the negative interactions that species in module III share with those in modules VII and VIII suggest that species in these groups limit the establishment of other species when they are already present in the patch.

Interestingly, the coarse-grained network of the Monegros ecosystem (Fig. 6.b) provides useful insights about the response of the system under certain perturbations. For instance, the aforementioned fragmentation into several groups of the population competing with group I can be seen as a product of the two unbalanced triads, (III-V-VIII) and (III-V-VII), in the coarse grained network. Thus, one can ask what are the effects associated to the removal of group III or, alternatively, the absence of species 33, as it is the one responsible of the unbalanced triads created in the network of species (Fig. 5.c). To this aim we have studied the Kuramoto dynamics with and without species 33 for K = 1. In this way, we evolve in time the Kuramoto dynamics of the whole system and at some point we remove node 33 in order to monitor how a new partition of the network shows up.

In Fig. 7 we report the original equilibrium distribution of the phases of the 44 species in the unit circle in the Monegros system (a) and the one after removing group III (b). The result is that the population passes from being partitioned into 10 groups to have 6 modules. One module (III) is trivially removed, however, as a byproduct other 3 groups have disappeared or, better said, three new groups have emerged from the mixing of six (IV-IX) old ones. The merging of the 6 modules is shown in Fig. 7.c by monitoring the evolution of ∆θ i for the nodes in groups IV-IX. At time t = 0 we remove group III and subsequently a reorganization of the species into new three groups takes place. The composition of the most populated of the new groups is highlighted in Fig. 5.c and, as shown, it comprises former groups IV, VI, VII and VIII and species 16 from group V. In a nutshell, the absence of species 33 creates a much larger coalition of species competing with those contained in group I. This result illustrates the critical role played by group III in the organization of the system and, in particular, in the fragmentation of the original coalition. In ecological terms, the suppression of species 33 would foster the formation of vegetation patches mixing species from former modules IV, VI, VII and VIII, that, in the original system, were highly conditioned by the presence of species 33.

VI. CONCLUSIONS

The analysis carried out in both Cabo de Gata and Monegros locations unveils that the interplay between negative and positive interactions yields an effective partition of signed net-works that can be detected by analyzing the synchronization patterns between the nodes. Importantly, this analysis can be done straightforward without any a priori information about the number and size of groups in the system and avoiding computationally expensive optimization algorithms.

The method introduced here is along the same line as the one reported in [START_REF] Arenas | Synchronization reveals topological scales in complex networks[END_REF][START_REF] Gómez-Gardeñes | Synchronizability determined by coupling strengths and topology on Complex Networks[END_REF][START_REF] Boccaletti | Detecting complex network modularity by dynamical clustering[END_REF][START_REF] Almendral | Dynamics of overlapping structures in modular networks[END_REF][START_REF] Sendiña-Nadal | Unveiling Protein Functions through the Dynamics of the Interaction Network[END_REF] for the case of community detection in complex (unsigned) networks by means of the classical attractive Kuramoto model. Very recently [START_REF] Wu | Partition signed social networks via clustering dynamics[END_REF], another dynamical clustering based on the dynamics of a deformed KM has been implemented to analyze the structure of signed social interactions. Here, after illustrating in Sec. IV the patterns that the interplay between attractive and repulsive interactions introduces in the stable equilibrium of Kuramoto oscillators in ideal modular signed networks [START_REF] Anderson | Multiscale dynamics in communities of phase oscillators[END_REF] by analyzing the case of two interconnected groups, in Secs. V.A and V.B we have tackled the problem of partitioning ecological signed networks with both competition and facilitation by simply analyzing the equilibrium distributions and compositions of the synchronization clusters. The identification of modules within the community and the position of species within those modules is central to understand the structure and functioning of ecological networks, independently of the type of biotic interactions considered [START_REF] Ings | Review: Ecological networks beyond food webs[END_REF].

The application of synchronization to group partitioning in ecological networks also provides useful insights about the nodes and groups which play a major role structuring the community. As shown in Sec. V.C, this information help us to explain the diversity and composition of vegetal patches and, more importantly, to anticipate the possible changes that the removal of species introduce in this organization. In this regard, the analysis of synchronization patterns allows to monitor how is the transition between the original organization to the emerging new one, thus providing a fast and useful framework to assess the resilience of ecological networks and to design accurate and focused conservation/restoration strategies.

FIG. 1 :

 1 FIG. 1: Two interconnected networks of 5 nodes each. Each network is a fully connected graph. The links among nodes of the same networks are positive (solid) whereas the 4 links connecting nodes of different networks are negative (dashed).

  FIG. 2: Panels (a)-(c) show the distribution of the values e iθ α i in the unit circle of the 2N (N = 500) oscillators for K = 0 (a), K = 0.5 (b) and K = 1 (c). The color of each circle denotes the natural frequency ω α i of each oscillator: the darker the more different from the average frequency (Ω = 0). In panel (d) we show the evolution of the synchronization of the system, r, and each of the networks, r α (α = 1, 2), as a function of the coupling K. The system is composed of two fully-connected networks interconnected with probability p = 0.4.

FIG. 3 :

 3 FIG. 3: Degree of synchronization within a population α as a function of K for different values of interconnection probability p. The value of K is normalized to that of the original Kuramoto model K KM c

A c c e p t e d m a n u s c r i p t FIG. 4 :

 t4 FIG. 4: Analysis of the Cabo de Gata ecosystem. The first panel (a) shows the evolution of the degree of synchronization r of the whole network as a function of the coupling strength K. Panel (b) shows the evolution, as a function of K, of the phase difference of node i with respect to node 1, ∆θi. Finally, panel (c) shows the network of interactions [facilitation links in black while negative ones appear in grey (orange in the color version)] by grouping the nodes in six groups as revealed from panel (b).

A c c e

  FIG. 5: Analysis of the Monegros ecosystem. A sin Fig. 4 the first panel (a) shows the evolution of the degree of synchronization r of the whole network as a function K. Panel (b) shows the evolution of ∆θi as a function of K, pointing out the existence of 10 groups. Finally, panel (c) shows the network of interactions [facilitation links in black while negative ones appear in grey (orange in the color version)] by grouping the nodes in the 10 groups derived from panel (b).

A c c e

  FIG.6: Coarse-grained description of the Cabo de Gata (a) and Monegros (b) ecosystems as obtained from the Kuramoto partitions. The nodes of these networks represent groups of species whereas a positive (negative) link between a pair of groups exists if, at least, one positive (negative) interaction between two species belonging to each group is present in the original network. Solid (dashed) links reprsent positive (negative) interactions.

A c c e p t e d m a n u s c r i p t FIG. 7 :

 t7 FIG. 7: Panels (a) and (b) show the distribution at equilibrium of the phases of each species in the unit circle. The panels show these distribution in the Monegros ecosystem before (a) and after (b) removing group III. Finally, panel (c) shows the time evolution of the values ∆θi of the nodes in groups IV-IX before (t < 0) and after (t > 0) removing group III.

A c c e

  p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

TABLE I :

 I Main properties of the networks analyzed. For each net-

	Area	N	⟨k⟩	⟨k⟩ +	⟨k⟩ -Nc S
	Cabo de Gata (Almería) 42 2.476 1.381 1.095 2 40
	Monegros (Aragón)	46 3.347 2.695 0.652 2 44

work (named as the area where data were collected) we show the number of nodes (species), N , the average degree of nodes, ⟨k⟩, the average number of positive (negative) connections per node, ⟨k⟩ + (⟨k⟩ -), the number of connected components, Nc, and the size of the largest one, S.

Acknowledgments

We are grateful to Concepción L. Alados and Yolanda Pueyo who kindly shared their data about Monegros (CGL2012-37508) and Cabo de Gata (CGL2011-27259). JGG acknowledges financial support from the European Commission through FET IP projects MULTIPLEX (Grant No. 317532) and PLEXMATH (Grant No. 317614), the Spanish MINECO through the Ramón y Cajal program and Projects FIS2015-71582-C2 and FIS2014-55867-P, from the Departamento de Industria e Innovación del Gobierno de Aragón y Fondo Social Europeo (FENOL group E-19) , and the Brazilian CNPq through the grant PVE of the Ciencias Sem Fronteiras program. HS acknowledges funding via a JAEPre grant from the National Research Council of Spain. RFSA and FSB acknowledge financial support from FAPESB (contract number PNX 0006/2009) and by CNPQ (INCT-SC).

APPENDIX A: LINEAR STABILITY OF THE SYNCHRONIZED SOLUTIONS IN THE TWO POPULATIONS SYSTEM

From the Kuramoto equations, Eq. ( 6), we can write the linearized equations governing the time evolution of small perturbations ϵ α i ≪ 1 to the phases θ α i as:

where δ ij is the Kronecker's delta function defined as δ ij = 1 if i = j and δ ij = 0 otherwise. Now we consider an initial perturbation of the form ϵ α i (0) = δ and ϵ β i (0) = -δ (with α = 1, 2 and β = 2, 1 respectively), i.e., we perturb all the phases in one population in the same way while the opposite perturbation is applied to all the phases in the other one. This way, the perturbations in population α initially evolve as:

whereas those in population β initially follow:

Let us note that we have used these two equalities:

Now we consider the first of the solutions to Eq. ( 8), Ψ α = Ψ β = 0 (θ α i = θ β i for those oscillators with the same identical frequency ω α i = ω β i ), and set (by invoking the symmetry of the two populations system) that r α = r β . Then we obtain from Eqs. (A2) and (A3): εα i = -εβ i > 0, so that the phases of the oscillators θ α i and θ β i separate from each other, thus pointing out that Ψ α = Ψ β = 0 is an unstable equilibrium.

Instead, substituting in Eqs. (A2) and (A3) the solution Ψ α = 0 and Ψ β = π (θ α i = θ β i -π for those oscillators with the same identical frequency ω α i = ω β i ), yields: εα i = -εi β < 0, so that the phases in both populations tend to recover their original position, i.e., they are stable.


Gautier Verhille
email: gautier.verhille@irphe.univ-mrs.fr

Adrien Bartoli

3D conformation of a flexible fiber in a turbulent flow

Keywords: Turbulent flows, Fibers in turbulence, Epipolar geometry PACS 47.27.-i, 42.30.Tz

)). We show that this function depends on an elastic length e which characterizes the particle flexibility, as is the case for polymers, but also on the fiber length L, contrary to polymers.

Introduction

The transport of fibers in turbulent flows is a key point in several applications such as the paper and textile industries [START_REF] Lundell | Fluid mechanics of papermaking[END_REF][START_REF] Marheineke | Fiber dynamics in turbulent flows: General modeling framework[END_REF] and polymer processing [START_REF] Jarecki | Modeling of pneumatic melt spinning processes[END_REF]. In all these cases, the quality of the final product depends on the fiber orientation statistics (homogeneity and isotropy mainly). The dynamics of fibers in turbulent flows is also relevant in environmental sciences. Indeed phytoplanktons play a key role in the ocean by producing organic matter and oxygen from the CO 2 dissolved in water. It is known that the rotation rate influences the chemotaxis of these organisms [START_REF] Macnab | The gradient sensing mechanism in bacterial chemotaxis[END_REF] and that the shape of a particle influences its rotation rate in a shear flow [START_REF] Hinch | Rotation of small non-axisymmetric particles in a simple shear flow[END_REF].

Since 2010, an increasing body of work has been devoted to the dynamics of anisotropic particles in turbulence. In theoretical works on anisotropic particles, the particle inertia is generally neglected. Then, the particles are assumed to have the same density as the carrying fluid and to be smaller than the Kolmogorov scale η K , i.e. the smallest spatial scale of the flow. In that case, the particles can be modeled as point particles with a given aspect ratio Λ (Λ ∼ ∞ for thin rod and Λ = 0 for thin discs). Thus the advection and the orientation dynamics are decoupled. While the advection is simply the one of a tracer following flow streamlines, the rotation dynamics is determined by the local shear according to Jeffery's work [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF]. As a consequence, fibers are preferentially aligned with the vorticity [START_REF] Pumir | Orientation statistics of small particles in turbulence[END_REF][START_REF] Chevillard | Orientation dynamics of small triaxial-ellipsoidal particles in isotropic turbulence[END_REF]. The 'point-fiber' approach has also been used in various numerical works such as in inhomogeneous and anisotropic turbulence to investigate the influence of walls on the fiber dynamics [START_REF] Marchioli | Rotation statistics of fibers in wall shear turbulence[END_REF]. The first measurement of the orientation dynamics of cylindrical particles in a turbulent flow was made by Parsa et al. [START_REF] Parsa | Rotation rate of rods in turbulent fluid flow[END_REF]. They used 4 high speed cameras to reconstruct the particle orientation in space. Another technique was used by Byron et al. [START_REF] Byron | Shape-dependence of particle rotation in isotropic turbulence[END_REF] with hydrogel particles seeded with small fluorescent beads. They could determine the global rotation rate of discs and rods with a classical Particle Image Velocimetry (PIV) algorithm. For more complex 3D objects, Marcus et al. [START_REF] Marcus | Measurement of the solid body rotation of anisotropic particles in 3d turbulence[END_REF] proposed to determine the orientation with a stereo-matching method, by minimizing the error between the projection of synthetic particles and experimental images. Zimmermann et al. [START_REF] Zimmermann | Tracking the dynamics of translation and absolute orientation of a sphere in a turbulent flow[END_REF] also used this technique to measure the rotation rate of a large painted spherical particle. The recorded pattern was then compared to a previous set of images to determine the particle orientation.

In all these works, the particles are assumed rigid. However, it is known that increasing the aspect ratio of a particle increases its flexibility. This raises the question of the influence of the deformation on the transport. In a previous paper, we showed that the flexibility of a fiber cannot be neglected if its length is longer than an elastic length e = (EI) 1/4 /(ρη ) 1/8 , where EI is the bending modulus of the fiber, ρ and η the fluid density and viscosity, and the power injected in the turbulence [START_REF] Brouzet | Flexibile fiber in a turbulent flow: a macroscopic polymer[END_REF]. This length is defined by the equilibrium between the injected power in turbulence and the power needed to bend a fiber. We also proposed an analogy between a fiber distorted by the turbulent fluctuations and a semiflexible polymer in thermal equilibrium in a good solvent. Measuring the 3D conformation of a fiber is a key point to test this analogy. But it may also be useful to determine some flow properties which can be difficult to estimate otherwise. For instance, the injected power in the turbulence can be known from the measurement of the elastic length e .

We present an experimental technique from 3D computer vision to reconstruct the 3D conformation of a fiber. By carefully adapting the method to our images we are able to determine the conformation of flexible fibers in a turbulent flow, from the images of two cameras. In this paper, we restrict ourselves to fiber lengths smaller or of the order of the flow integral scale and to a single elastic length which lies in the inertial range. We focus on the evolution of the correlation function of the orientation of the tangent vector along the fiber, which is a key quantity to understand the analogy with polymers. We will see that the shape of this function is similar in both cases, but that the particle length plays a role for fibers contrary to polymers. The experimental setup is described in the second section. The third section is focused on the 3D reconstruction technique and the fourth section presents the validation of the technique. Finally, our results are presented and discussed in the last section, before the conclusions.

Experimental setup

Turbulence is generated in a von Kármán flow by the counter-rotation of two impellers [START_REF] Ravelet | Bifurcations globales hydrodynamiques et magnétohydrodynamiques dans un écoulement de von Kármán[END_REF][START_REF] Zimmermann | Tracking the dynamics of translation and absolute orientation of a sphere in a turbulent flow[END_REF]. Whereas most experiments are done in a cylindrical geometry, the flow is here confined in a cubic box of 20 × 20 × 20 cm 3 to reduce optical distorsions, cf. a diameter of 2R=17 cm and are fitted with 6 straight blades with height 0.5 cm. The distance between the two impellers is 17 cm. They are set in motion by two 350 W asynchronous motors at a constant rotation frequency F in the range [1;30] Hz. All the results presented here have been measured in water at ambient temperature (ρ = 10 3 kg.m -3 , η = 10 -3 Pa.s) with a rotation frequency of 15 Hz. The energy dissipation rate is then around 6 W.kg -1 . This flow has been extensively used to study different problems in turbulence such as Lagrangian turbulence (see for instance [START_REF] Toschi | Lagrangian properties of particles in turbulence[END_REF] and references therein), dynamo instability [START_REF] Miralles | Dynamo threshold detection in the von Kármán sodium experiment[END_REF] or quantum turbulence [START_REF] Rousset | Superfluid High REynolds von Kármán experiment[END_REF]. The mean flow is composed of two toroidal cells with a large azimuthal shear at the center and two poloidal recirculation cells due to the centrifugal forces near the impellers [START_REF] Machicoane | Large sphere motion in a nonhomogeneous turbulent flow[END_REF][START_REF] Lopez-Caballero | Large scales in a von Kármán swirling flow[END_REF]. If the small scales of the flow are expected to be isotropic due to the high Reynolds number, the large scales are anisotropic due to the stretching induced by the poloidal cells [START_REF] Voth | Measurement of particle accelerations in fully developed turbulence[END_REF][START_REF] Ouellette | Smallscale anisotropy in lagrangian turbulence[END_REF]. The influence of the flow anisotropy is beyond the scope of this study and will be addressed in future work.It is known that the turbulence homogeneity is fairly good close to the stagnation point and decreases as one approaches the impellers. That is why the reconstruction volume is a cube with approximately 6 cm side length centered at the stagnation point of the flow (see the black cube in figure 1).

As discussed in the introduction, the fiber deformations are important when the fiber length L is larger than the elastic length e . Due to the size of the reconstruction volume we restrict ourselves to fiber length smaller than L max = 5 cm. Then, to maximize the deformation, we have to reduce the elastic length e = (EI) 1/4 /(ρη ) 1/8 . Working with water at a rotation frequency of 15 Hz determines the denominator. Then the only free parameter to minimize this characteristic length is the bending modulus, i.e. the Young modulus E and the fiber diameter d (the area moment of inertia I is proportional to d 4 ). Softer materials with very low Young modulus such as hydrogel are generally fragile. They cannot be used to cast fibers as they would be destroyed when hitting an impeller. We found a good compromise with the silicone EC00, from Esprit Composite, which has a Young modulus E = 40 ± 15 kPa. Even if this incertitude is relatively large it does not affect the elastic length significantly as it scales as E 1/4 . Then the uncertainty on the elastic length is of the order of 10 %. Since this silicone is very extensible (more than 400 % of deformation before breaking), fibers are not broken when an impact occurs. This translucent silicone has been dyed with a manganese colorant powder (particle size between 1-2 µm) to improve the contrast for the particle detection. To make fibers the silicone is first molded in a straight laiton tube with a diameter of 600 µm and then cut to the desired length, here between 1 and 5 cm. We have checked the fiber diameter on several fibers with a microscope and found d = 622 ± 13 µm. With these properties, the estimated elastic length e is around 3.4±0.3 mm which is nearly 1/3 of the length of the shortest fibers. One can then expect that all fibers will be greatly distorted by the flow. As the silicone density is 1.03, the fiber settling speed is of the order of 1 cm.s -1 . This velocity has to be compared to the flow speed which is of the order of 1 m.s -1 . Therefore buoyancy effects can be neglected in the problem. Finally, to avoid inter-particle interactions the total volume concentration of fibers is around 0.002 % (so between 3 and 20 fibers depending on their length).

Since this paper focuses on the statistics of the fiber deformations and not on its dynamics, the use of high speed cameras is not necessary. Furthermore, to converge the statistics with less data, the image acquisition rate is low (a few images per second) to record uncorrelated conformations. However, the presented technique could be applied directly to dynamical problems such as the relaxation of the elastic strain or the existence of elastic waves propagating along the fiber. Two IDS UI-5240CP 1.3 Mpixels cameras with 20 mm lens are used to film the fibers in the flow at two different angles. The images are then processed with Matlab to determine the 3D shape of the fiber. To increase the contrast we used back-lighting with a LED-panel and a LED spot put behind diffusers. These light sources allow us to image fibers with an exposure time of 0.3 ms, chosen to avoid motion blur, and at a frequency of 5 images per second. This ensures that images are uncorrelated as the acquisition rate is 3 times smaller than the rotation frequency of the impellers. The angle between the two axes of the cameras is nearly 90 • and is determined precisely during the calibration process described in section 3.1.

3D reconstruction of the fiber conformation

Principle

We intuitively know that given two images, one can reconstruct the observed 3D structure. This intuition is based on the principle of triangulation, whereby the sightlines associated to two corresponding points, one in each image, meet in the 3D space where the sought 3D point lies. This intuition forms the backbone of 3D computer vision, and has been formalized and thoroughly studied in the past few decades [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Faugeras | The Geometry of Multiple Images[END_REF]. 3D computer vision has numerous applications, such as depth perception in computer-aided surgery [START_REF] Maier-Hein | Comparative validation of single-shot optical techniques for laparoscopic 3D surface reconstruction[END_REF].

Modeling one camera. The most popular camera model used in 3D computer vision is the so-called pin-hole camera. It abstracts the sensor and lens forming the camera by a simple plane, onto which the image is formed by intersecting a line joining the 3D point and the center of projection of the camera. Optical distortions, especially radial and tangential distortion, are then modeled directly in the image space. In the pinhole model, a camera has 6 extrinsic and 5 intrinsic parameters. The extrinsic parameters give the position and orientation of the camera in space and are represented by a 3D rotation R and translation t in some world coordinate frame. The intrinsic parameters give the focal length f of the camera, the principal point (the orthogonal projection of the center of projection onto the image plane) with coordinates c x , c y and the skewness τ and aspect ratio α of the pixels. They are held in an upper triangulation 3 × 3 matrix K whose bottom-right entry is fixed to one:

K =   αf τ c x 0 f c y 0 0 1   . (1)
As will be seen shortly, this representation allows one to easily project 3D points to the image using homogeneous coordinates. More precisely, let i ∈ {1, 2} be the camera index, we have the extrinsic parameters as R i , t i and the intrinsic parameters as K i . We now define the 3 × 4 projection matrix P i as:

P i ∝ K i   1 0 0 0 0 1 0 0 0 0 1 0   S i , with S i =     R i t i 0 0 0 1     . ( 2
)
The proportionality arises from the use of homogeneous coordinates. In that system of coordinates two points

Q 1 (X 1 , Y 1 , Z 1 , H 1 ) and Q 2 (X 2 , Y 2 , Z 2 , H 2 ) are identical (Q 1 ∝ Q 2 ) if: (X 1 , Y 1 , Z 1 , H 1 ) = λ(X 2 , Y 2 , Z 2 , H 2 ), λ ∈ R . (3)
For a 3D point whose world coordinates are (X, Y, Z) , we construct its homogeneous coordinates, defined up to scale, as (X, Y, Z, 1) . The advantage of these coordinates is that the projection operation is then simply defined as q i ∝ P i Q, where q i represents the homogeneous coordinates of the projected point in image i.

Modeling two cameras. The geometry of two cameras is specifically called the epipolar geometry [START_REF] Longuet-Higgins | A computer program for reconstructing a scene from two projections[END_REF]. For the pin-hole camera model, the epipolar geometry is represented by the 3 × 3 fundamental matrix F . This matrix can be constructed from the two projection matrices P 1 , P 2 . Given a point with homogeneous coordinates q 1 in the first image, the fundamental matrix associates the homogeneous coordinates of a line l 2 ∝ F q 1 to it in the second image. This line is called the epipolar line, and represents the projection of the sightline of q 1 for the first camera into the second camera. It is obvious that the corresponding point q 2 of q 1 in the second image must lie on l 2 . We thus have the following algebraic constraint, called the fundamental equation: q 2 • F q 1 = 0. We will use this constraint in the next section to match the two images of a fiber, finding corresponding points by intersecting the fiber's detected image curve with epipolar lines. The epipolar lines can obviously be constructed in the first image from the second image, and this is done by using

l 1 ∝ F q 2 .
Calibrating a pair of cameras. The process of calibration determines the intrinsic and extrinsic parameters of the two cameras, from which the two projection matrices P 1 and P 2 are formed, as well as the fundamental matrix F . Calibration is done by showing the camera an object with known structure. In practice, the calibration object is constructed by printing a checkerboard pattern using a standard laser printer, which is then plasticized to make it waterproof and to increase its rigidity. Several image pairs of the calibration object are recorded with various orientations, and calibration is solved by Bouguet's method [START_REF] Bouguet | Camera calibration toolbox for Matlab[END_REF]. This process also estimates the radial and tangential distortion parameters, which can be neglected in our case. These will be systematically used to 'undistort' any image collected by the cameras from then on. In practice, it is important that the image pairs of the calibration object are recorded in situ. Indeed, this will guarantee that the calibration is valid where it needs to be, in other words, where we will need to reconstruct the fibers' 3D conformation. It is also important because our experimental setup has two interfaces, between air and Plexiglas, and Plexiglas and water. This raises the question of adapting the pin-hole and epipolar models to this specific configuration, as was studied in [START_REF] Chari | Multi-view geometry of the refractive plane[END_REF]. It was however showed in a recent paper that calibrating in situ naturally adapts the models' parameters to cope with this case very well, especially for medium field-of-view cameras [START_REF] Treibitz | Flat refractive geometry[END_REF].

Triangulating points. We are now back to our initial intuition that 3D points can be reconstructed by intersecting two corresponding sightlines from cameras with distinct centers of projection. Concretely, triangulation is the process of finding Q from q 1 , q 2 , P 1 , P 2 .

The camera matrices P 1 , P 2 were found by calibration, and the corresponding points q 1 , q 2 are found by the matching process described in the next section. Algebraically, the triangulation constraints are written as q 1 ∝ P 1 Q and q 2 ∝ P 2 Q. One notices that each projection gives two constraints (as the predicted image point P i Q must match the two coordinates of the measured image point q i ), and so a total of four constraints for the two images, but that we have only three unknowns in the 3D point Q. Geometrically, this is because one has to intersect the two sightlines, which may not meet. The most successful approach is to find the 3D point whose image projections are closest to the measured image points. This problem has an analytical solution which was given in [START_REF] Hartley | [END_REF].

Fiber reconstruction

This section is devoted to the description of how the 3D fiber conformation is determined from two images. In elasticity, the deformation of a slender object is defined by its neutral line. For an homogeneous cylinder, as is the case here, the neutral line corresponds to the center of the cross section. To determine the position of this line, two approaches are possible. The first one consists in reconstructing directly the fiber centerline by an estimation in each image of its projection, and a triangulation of these points. In the second method, the fiber edges are determined in each image and then triangulated in space. The neutral axis is then determined from the position of the edge points in space. Both techniques have been tested and give similar results for weakly distorted fibers. However if one projection of a fiber intersects itself or presents a large curvature, the approach based on the centerline frequently gives inaccurate results preventing the fiber reconstruction. For these cases the edge technique seems more robust.

We will thus only discuss the results obtained with this method in the following. First, for each image, the edge coordinates are determined to sub-pixel precision using a Matlab function developed by Trujilo et al. [START_REF] Trujilo-Pino | Accurate subpixel edge location based on partial area effect[END_REF]. This routine returns the point coordinates and the local normal vector to the edge. Each fiber is then individualized by connecting each point to its nearest neighbors. If two fibers intersect on an image they are considered as a single object after this step. Then, for each object, the fiber extremities are determined using the approach given in Lee [START_REF] Lee | Curve reconstruction from unorganized points[END_REF] and matched from one image to the other. The edges are then oriented starting from the matched extremities in order to define a curvilinear coordinate system. To define the orientation one starts from the extremity and moves from closest neighbor to closest neighbor. This is not possible if there is an overlap of two parts of the same fiber or of two different fibers. To overcome this, we calculate the angle between the mean tangent vector calculated over the last 10 points and the vector linking the current point to its closest neighbor. If this angle is larger than some maximal value, we look for the closest point which is in the direction of the mean tangent vector and with the local normal vector aligned with the mean local normal vector calculated over the last 10 points. An illustration of the sorting of the edge points is presented in figure 2. Then each point of the edge of the fiber in image 1 is matched with the fiber edge points in image 2. Since the fiber is distorted, the epipolar line may intersect the contour at several locations. However, since the contour is oriented it is possible to eliminate false matched points by moving from closest neighbor to closest neighbor in image 1. In that case the matching point is the closest point to the previous point in image 2.

Once we found the two points q 1 and q 2 on each image, we determine the 3D coordinates of point Q by triangulation. At the end of this process, we obtain a set of 3D data points representing the fiber. We then need to define the centerline coordinates. First the curvilinear coordinate s is defined for each points Q by the length of the shortest path between Q and the matched extremity used previously. This shortest path is determined by making a graph from the point coordinates defining the weight of each link by the distance between points. Links larger than the fiber diameter are removed in order to isolate each fiber. Then the centerline coordinate is found by fitting the data points with a cubic spline function. The smoothing parameter is chosen so that the maximum curvature of the reconstructed fiber does not exceed the inverse of the fiber diameter, as is the case in linear elasticity.

Validation

In practice a fiber might not be entirely in the volume seen by both cameras. In that case, only a part of the fiber can be reconstructed. It is not possible to use this fiber for the statistics as the fit imposes that the fiber curvature is equal to zero at the fiber extremities. To eliminate these data, we only keep fibers with an estimated length L e of the order of the real fiber length L, in practice L e = L ± 10 %.

An illustration of the reconstruction is represented in figure 3 where the plain line is for the fitted curve and the scattered data for the raw data points. Note that for clarity, the distance between the raw data points and the fitted curve has been multiplied by 2. The validity of the reconstruction can be first emphasized by comparing the estimated length L e with the fiber length L. It has to be noticed that the way fibers are manufactured gives an uncertainty on L of the order of 1 mm. One Fig. 3 Exemple of fiber reconstruction: gray points correspond to raw data from the 3D reconstruction and the red curve is the cubic spline fit. Exemple of fiber reconstruction: the gray points correspond to raw data from the 3D reconstruction and the red curve is the cubic spline fit. For clarity, nearly 70 % of raw data points were removed randomly and the distance between the centerline and the data was multiplied by 2.

eraging depends on the fiber length. There are two reasons for that. First, the fiber concentration was smaller for long fibers to facilitate the fiber reconstruction and prevent the fiber interaction. Then, longer fibers were trapped close to the impellers more frequently, as observed for a large rigid sphere by Machicoane et al. [START_REF] Machicoane | Large sphere motions in a nonhomogeneous turbulent flow[END_REF]. Therefore, longer fibers enter in the reconstruction volume less frequently. However, we have checked that the amount of conformation was enough to converge the different statistics we present in each cases. In table 1, we can also see that the estimation of the mean values of the estimated length L e , where • stands for an ensemble average and the fiber length are equal within the manufactured precision. This result may be trivial as only fibers in the range [0.9 ; 1.1]L are selected for the analysis. However the distribution of the measured length is strongly peaked. This is shown by the value of the ratio of the standard deviation over the mean of the estimated length L rms e / L e which is nearly twice as small as the one of a randomly distributed variable. This also shows that, as for laminar flows, the stretching of the fiber is negligible as compared to the bending [START_REF] Powers | Dynamics of filaments and membranes in a viscous fluid[END_REF] and that the fiber can be considered as inextensible.

In addition, the technique used here allows us to determine the fiber diameter. Indeed, as the two camera axes are perpendicular, the edges of the projec-tion of one image correspond to the midline of the fiber projection in the second image. Then the triangulation of the intersection of the epipolar line with the fiber edge defines a square with side length equal to d, see figure 4. Then the fiber diameter is mea- sured by determining the distance between the reconstructed points and the fitted line. We found here d m = 626 ± 30 µm. As mentioned in section 2, the fiber diameter has been measured independently with a microscope to d = 622 ± 10 µm. These two measurements are in very good agreement.

In this section we have shown that the reconstructed points are at a distance d/ √ 2 of the neutral line fitted with the spline function. We have also seen that the length of the fitted spline is close to the real fiber length. These observations validate the reconstruction process. We now focus on the evolution of the tangent vector along the fiber.

Results and discussion

In Brouzet et al. [START_REF] Brouzet | Flexibile fiber in a turbulent flow: a macroscopic polymer[END_REF], an analogy between a wormlike chain polymer and a flexible fiber in a turbulent flow has been drawn. Indeed the evolution of the norm of the end to end vector R e , linking the two extremities, as a function of the fiber or polymer length seems similar in both cases. In the wormlike chain model [START_REF] Yamakawa | Modern Theory of Polymer Solutions[END_REF] this norm depends on two independent parameters, the polymer length L and the persistence length l p :

R 2 e = 2Ll p -2l 2 p 1 -e -L/lp . (4)
The persistence length l p is related to the polymer rigidity and is defined by EI/l p = kT where EI is the bending modulus and kT the thermal energy.

In fact in polymer theory, this relation is an integral relation and is derived from the correlation function of the orientation of the tangent vector:

C( ) = t(s) • t(s + ) s = e -/lp . ( 5
)
Here t is the tangent vector (with t • t = 1), s and are 2 curvilinear coordinates and • s represents a double averaging process over realizations and over the curvilinear coordinate s. We propose to measure this correlation function for different fiber lengths in the turbulent von Kármán flow described earlier. The results are then compared to the wormlike chain prediction. We first focus on the correlation function for the longest fibers (L = 5 cm). This function is presented in figure 5 where the distance has been normalized by the elastic length scale e . In this graph the experimental curve (plain blue line) is compared to the exponential law predicted for wormlike chain polymers, cf. equation [START_REF] Chevillard | Orientation dynamics of small triaxial-ellipsoidal particles in isotropic turbulence[END_REF]. The global evolution is well captured by this prediction with a persistence length l p 7.5 e . However the inset of figure 5 shows that at short length the evolution of the experimental curve is quadratic (C( ) ∼ 1 -α 2 ) whereas for an exponential law, the evolution is linear (exp(-/l p ) ∼ 1 -/l p ).

In fact, one can easily show that this linear behavior is not possible either for fibers or for wormlike chains. Indeed the Taylor expansion of the correlation function gives:

t(s)•t(s + ) t(s) • t(s) + t(s) • ∂ s t(s) + 2 2 t(s) • ∂ ss t(s) + 3 6 t(s) • ∂ sss t(s) + O( 4 ). (6)
By definition the tangent vector is normalized. Thus it has to verify the two relations t • t = 1 and ∂ s (t • t) = 0. Therefore for small , the correlation of the orientation can be written as:

t(s)•t(s + ) 1 -κ 2 2 /2 -∂ s κ 2 3 /4 + ..., (7)
in contradiction with the simple exponential approximation. Here κ 2 = ∂ s t • ∂ s t is the local curvature of the fiber. As this relation does not involve any physical properties of the fiber the same behavior is expected for polymer. Furthermore this relation enlightens the fact that the persistence length is a function of the local curvature of the fiber. This idea is already present in Brouzet et al. [START_REF] Brouzet | Flexibile fiber in a turbulent flow: a macroscopic polymer[END_REF]. They explained that if the fiber length L is equal to the elastic length e , the correlation function should be equal to sin θ/θ, where θ = / e . Indeed, when the fiber length is equal to the elastic length, the deformation is an arc of circle with radius R e . Then, the correlation function is simply:

t(s) • t(s + ) = cos /R. (8)
The correlation function can be found by averaging over all configurations, i.e. from a straight line to an arc of circle of radius e :

C( ) = t(s) • t(s + ) = L sin L = e sin e . (9)
We have compared this prediction to our measurements in the inset of figure 5. We can see that the cardinal sine (dashed dotted line) fits very well the experimental curve with no adjustable parameters for < 0.5 e . Therefore, the Taylor expansion of the cardinal sine (sin(x)/x ∼ 1 -x 2 /6) gives a relation between the elastic length e and the mean curvature κ 2 : We focus now on the evolution of the correlation function with the fiber length for the same elastic length e . We can see in figure 6 that the experimental curve does not superimpose for all lengths. Whereas the correlation functions overlap for the shortest and the longest fibers, they decrease more rapidly for intermediate lengths. This is in contradiction with the wormlike chain model which assumes that the correlation function depends only on the persistence length l p . However, for long enough fibers (L > 1 cm) two regimes are visible. For large , the correlation function is nearly linear in a log-lin plot (see figure 6). So this part can still be modeled by an exponential law but with a characteristic length which depends also on the fiber length. For small , the correlation function is still quadratic but the mean curvature κ = ( κ 2 ) 1/2 changes with the fiber length. Then the fit with the cardinal sine seems to be valid only for the shortest and the longest fibers. The relation between the mean curvature and the elastic length does not hold for the intermediate lengths. An evaluation of κ with the elastic length using equation (10) will underestimate the value by nearly 10 % for the 2 cm fiber.

κ 2 = 1 3 2 e . (10)
As the shape of the correlation function varies with the fiber length in turbulence, the analogy between the deformations of a fiber in a turbulent flow and the distorsions of a polymer submitted to thermal fluctuation is not rigorous. We may now wonder if this result could also have been observed without the 3D reconstruction tools we have developed. In Brouzet et al., the analogy was based on the non linear evolution of the norm of the end to end vector R 2 e when the fiber length is increased. The same quantity is presented in figure 7 where the measurement uncertainty is nearly the symbol size. The dashed line represents the best fit of the evolution of the norm of the end-to-end vector as a function of the fiber length, cf. equation ( 4), where the persistence l p is considered here as a fit parameter and is noted l p,Re in the following. We can see here that with only this integral quantity it is difficult to assess the analogy between polymers and flexible fibers as the fit with l p,Re 6 e is very close to the experimental points. The 3D reconstruction allows one to access local quantities, such as the correlation function, which are more sensitive and allowed us to clearly show that the fiber length plays an important role in the fiber deformation in turbulence contrary to wormlike chain polymer.

Note that in Brouzet et al. the persistence length was equal to 2 e and not 6 e as is the case here. This difference is not due to the difference of the experimental techniques used in both studies but is related to the flow correlation. Indeed, in turbulence the flow is correlated at scale r smaller than the so called integral scale L I with a correlation C t (r) ∼ r 2/3 . As the elastic length is here 10 times smaller than in our previous study, the turbulent fluctuations at the origin of the deformation do not have the same correlation. Then it is expected that the relation between the persistence length and the elastic length differs in both studies. Beyond the polymer analogy, this observation opens new perspectives to model the fiber deformation in turbulence which appears to be a very rich problem.

Conclusion

In this study we used the classical epipolar geometry to determine fiber conformations in a turbulent von Kármán flow. The measurement of the correlation function of the orientation of the tangent vector along the fiber shows that the analogy between flexible fibers distorted by turbulent fluctuations and wormlike chain polymers submitted to thermal fluctuations is not always true. In both cases, the correlation function C(s) is proportional to 1-κ 2 s/2 for s e and, for s e , it can be approximated by an exponential law exp(-s/l p ). However we showed that for flexible fibers, the persistence length is a function of both the elastic length e and the fiber length L whereas it depends only on the elasticity for polymers. This dependence on the fiber length is directly related to the flow correlation, which is scaledependent in turbulence.

The relation between the flow properties and the fiber statistics deserves to be studied in more details. As is shown by equation [START_REF] Hartley | [END_REF], there is a direct relation between the correlation function and the local curvature. A better understanding on the physics of the deformation might be possible by investigating the curvature statistics and the tools developed here can be easily adapted to this task. Finally, the measurement of fiber deformations might allow one to determine some flow properties. For instance, the estimation of the persistence length for very long fibers seems related to the flow integral scale. It is also possible to estimate the turbulent dissipation rate by looking at the correlation function for short distances. Indeed, by fitting the correlation function with a quadratic function, it is possible to estimate the elastic length e = (EI) 1/4 /(ρη ) 1/8 within 10 %. This gives within a relatively large range as it can be overrestimated by a factor of 2. This kind of measurement could however be useful in oceanography as it allows a rapid estimation of the local dissipation rate.

figure 1 .Fig. 1

 11 Fig. 1 Schematic of the experimental setup. The black cube in the center of the vessel materializes the reconstruction volume where the fiber conformation is measured.

Fig. 2

 2 Fig.2Raw images of distorted fibers in a turbulent flow. The result of the contour detection and of its orientation is presented in the inset where the crosses represent the edge points and the color, from red to yellow, represents the curvilinear coordinate (we display only one half of the points for clarity).

  the plane containing both camera axis

Fig. 4

 4 Fig.4Triangulation of points from the edges of a cylinder. The red (resp. green) line represents the edges measured in image 1 (resp. 2). The triangulated points are at the intersection of the sightlines from cameras 1 and 2. As shown at the bottom of the image, they define a square tangent to the cylinder with side length equal to d.

Fig. 5

 5 Fig.5Evolution of the correlation function for the experimental data in plain (blue) line, the best exponential fit (black, dashed line). In the inset, zoom for small of the correlation function (plain line) and the dashed dotted line represents the cardinal sine function sin(x)/x, cf. text for details.

Fig. 6

 6 Fig.6Evolution of the correlation length as a function of the fiber length.

2 ]Fig. 7

 27 Fig.7Mean value of the norm of the end to end vector R 2 e as a function of the fiber length. The blue dot is the measurement and the dashed line the best fit from the prediction of the wormlike chain theory. The uncertainty scales roughly with the symbol size.

Table 1

 1 Comparison of the measured length and the fluctuation rate for all fiber length. The last row indicates the number of fibers used to compute the different presented statistics.

	Fiber length: L [cm]	1	2	3	5
	Estimated length: Le [cm]	1.01	1.92	2.90 5.16
	Fluctuation rate: L rms e	/ Le [%]	3.8	3.6	3.2	3.0
	Number of conformations	8868 1501	791	296

can see in table 1 that the number of fibers used for av-

Acknowledgements This work has been carried out in the framework of the Labex MEC Project (No. ANR-10-LABX-0092), of the A*MIDEX Project (No. ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir" French Government program managed by the French National Research Agency (ANR). A. Bartoli was funded by the FP7 ERC research grant 307483. Authors want to thanks B. Favier and P. Le Gal for fruitful conversation.