P-P Cortet

E Herbert

A Chiffaudel
email: arnaud.chiffaudel@cea.fr

F Daviaud

B Dubrulle

V Padilla

Susceptibility divergence, phase transition and multistability of a highly turbulent closed flow

Keywords: Turbulence, Metastable states, Critical exponents and amplitudes (Experiments), Classical phase transitions (Experiments)

Using time-series of stereoscopic particle image velocimetry data, we study the response of a turbulent von Kármán swirling flow to a continuous breaking of its forcing symmetry. Experiments are carried over a wide Reynolds number range, from laminar regime at Re = 10 2 to highly turbulent regime near Re = 10 6 . We show that the flow symmetry can be quantitatively characterized by two scalars, the global angular momentum I and the mixing layer altitude z s , which are shown to be statistically equivalent. Furthermore, we report that the flow response to small forcing dissymetry is linear, with a slope depending on the Reynolds number: this response coefficient increases non monotonically from small to large Reynolds number and presents a divergence at a critical Reynolds number Re c = 40 000 ± 5 000. This divergence coincides with a change in the statistical properties of the instantaneous flow symmetry I(t): its pdf changes from Gaussian to non-Gaussian with multiple maxima, revealing metastable non-symmetrical states. For symmetric forcing, a peak of fluctuations of I(t) is also observed at Re c : these fluctuations correspond to timeintermittencies between metastable states of the flow which, contrary to the very-longtime-averaged mean flow, spontaneously and dynamically break the system symmetry. We show that these observations can be interpreted in terms of divergence of the susceptibility to symmetry breaking, revealing the existence of a phase transition. An analogy with the ferromagnetic-paramagnetic transition in solid-state physics is presented and discussed.

Introduction

Symmetry breaking is an essential ingredient of classical phase transitions [START_REF] Landau | Statisticheskaya Fisika Third edition[END_REF]. Remarkably, it also governs the transition to turbulence, that usually proceeds, as the Reynolds number Re increases, through a sequence of bifurcations breaking successively the various symmetries allowed by the Navier-Stokes equations coupled to the boundary conditions [START_REF] Manneville | Dissipative Structures and Weak Turbulence[END_REF]. Finally, at large Reynolds number, when the fully developed turbulent regime is reached, all the broken symmetries are restored in a statistical sense [START_REF] Frisch | Turbulence -The Legacy of A N Kolmogorov[END_REF].

Beside this well-established scenario, some experiments in closed flows have raised intriguing features proceeding during transition to turbulence [START_REF] Chillá | [END_REF]5,6,7,8,9,10,11]. For instance, Tabeling et al. have evidenced a local peak in the flatness of velocity derivative in a von Kármán flow of helium around Re = 2 × 10 5 [8] that was claimed to exhibit some characteristics of a second-order phase transition [9] and suggested to be associated with the breakdown of small-scale vortical structures. However, since then, no numerics or experiments ever confirmed this scenario, leaving this flatness peak as an unsolved problem. More recently, a transition has been evidenced in another von Kármán flow around Re = 10 4 [10,11] which consists in a global bifurcation of the mean flow from the basic symmetric flow topology towards two other flow topologies which spontaneously break the symmetry of the driving apparatus. This transition is associated with an hysteresis and a divergence of transition times, two features classically observed in phase transitions. Even more recently, in a von Kármán flow very similar to the flow considered here, de la Torre and Burguete observed a time-intermittency of the large Reynolds number turbulent state, switching in between two metastable symmetrybreaking states [12,13]. Such transitions with dynamical symmetry breaking of mean patterns of the turbulent state are observed in 2D turbulence simulations [14,15], where they have been claimed to have strong geophysical and astrophysical relevance. They have also been observed and modeled in dynamo regimes of a liquid-sodium turbulent flow [16,17].

Turbulent flows being intrinsically out-of-equilibrium systems, there is a priori no reason to describe them using tools borrowed from thermodynamics. However, in recent years, developments of out-of-equilibrium thermodynamics led, in particular, to a number of interesting applications to turbulence [18,19,20]. In this context, one may wonder whether the observed turbulent transitions can be interpreted in terms of phase transitions with a symmetry-breaking or susceptibility divergence signature. In the above mentioned works, while the Reynolds number appears as the natural control parameter, the equivalent of an order parameter and its associated susceptibility has not been clearly identified.

In this paper, we consider a von Kármán turbulent flow and its averaged largescale patterns. We introduce a susceptibility to symmetry breaking for this flow and investigate its evolution as Re increases from 10 2 to 10 6 using Stereoscopic Particle Image Velocimetry (SPIV). We show that this susceptibility is as much as 50 times higher for the turbulent flow at very high Reynolds number than for the laminar flow. Furthermore, in a intermediate Reynolds number range -20 000 Re 200 000-, we measure susceptibilities values as high as 300 times the laminar value and conclude to the existence, at a critical Reynolds number Re c = 40 000 ± 5 000, of a divergence of this susceptibility to symmetry breaking. We also show that this phenomenon is associated with a peak in the amplitude of the fluctuations of the flow instantaneous symmetry. These observations complete the results of reference [21] with new complementary measurements at Re Re c and with a re-processing of all the time-series. These new statistical analysis allow us to consider not only the mean values of the distributions but also their most-probable values. In particular, the divergence of the susceptibility and the peak of fluctuations which were previously found to occur at different values of Re [21] are now observed at the same Reynolds number if we consider the susceptibility based on most-probable values instead of mean values. This result sets the threshold for a "turbulent phase transition" at critical Reynolds number Re c = 40 000 ± 5 000. We relate the especially large fluctuations observed close to the transition to timeintermittencies associated with an ergodic exploration of a band of metastable states with spontaneous symmetry breaking.

System, flow, symmetries and order parameters

Experimental setup

Our experimental setup consists of a Plexiglas cylinder of radius R = 100 mm filled up with either water or water-glycerol mixtures (cf. figure 1(a)). The fluid is mechanically stirred by a pair of coaxial impellers rotating in opposite sense. The impellers are flat disks of radius 0.925 R, fitted with 16 radial blades of height 0.2 R and curvature radius 0.4625 R (cf. figure 1(b)). The disks inner surfaces are 1.8 R apart setting the axial distance between impellers from blades to blades to 1.4 R. The impellers rotate, with the convex face of the blades pushing the fluid forward, driven by two independent brushless 1.8 kW motors. The rotation frequencies f 1 and f 2 can be varied independently from 1 to 12 Hz, speed servo loop control ensuring a precision of 2 for the mean frequency f = (f 1 + f 2 )/2 and, for the relative frequency difference θ = (f 1 -f 2 )/(f 1 + f 2 ), an absolute precision of 1 × 10 -3 and a stability along time better than 0.5 × 10 -3 .

Velocity measurements are performed with a SPIV system provided by DANTEC Dynamics. The cylinder is mounted inside a water or glycerol filled square Plexiglas container in order to reduce optical deformations. Two digital cameras are aiming at a meridian plane of the flow through two perpendicular faces of the square container providing 2D-maps of the three-components velocity field. Correlation calculations are performed on 32 × 32 pixels 2 windows with 50% overlap. As a result, each velocity is averaged on a 4.16 × 4.16 mm 2 window over the 1.5 mm laser sheet thickness. The spatial resolution is 2.08 mm. The data provide the radial u r , axial u z and azimuthal u ϕ velocity components on a 95×66 points grid through time series of 400 to 27 000 fields regularly sampled, at frequencies from 1 to 15 Hz, depending on the turbulence intensity and the related need for statistics.

The control parameters of the studied von Kármán flow are:

(i) the Reynolds number Re = 2πf R 2 /ν, where ν is the fluid viscosity, which controls the intensity of turbulence and,

(ii) the rotation number θ = (f 1 -f 2 )/(f 1 + f 2 )
, which controls the asymmetry of the forcing conditions.

System symmetry vs. flow symmetry

With exact counter-rotation of the impellers, i.e. when θ = 0, the experimental system (cf. figure 1(a)) is symmetric with respect to any R π -rotation exchanging the two impellers: the problem conditions are invariant under π-rotation around any radial axis passing through the center of the cylinder (cf. figure 1(c)). The symmetry group for such experimental system is O(2) [22]. When the motors frequencies differ, i.e. when θ = 0, the experimental system is no longer R π -symmetric and the symmetry switches to the SO(2) group of rotations. However, the parameter θ, when small but non-zero, can be considered as a measure of the distance to the exact O(2) symmetry: the stricto sensu SO(2) system at small θ can be considered as a slightly broken O(2) system. Such perturbative approach of the symmetry breaking has been successfully applied [23,24] for the 1:2 spatial resonance (or k -2k interaction mechanism) with slightly broken reflection symmetry. Driven at a given θ and a given Re, the system produces a flow. This flow -or at least its time-average at high Re-may be -or may not be-R π -symmetric. The problem is now to define a scalar quantity S, equivalent to an order parameter, which can quantify the distance of the flow to the R π -symmetry. The evolution of S with θ will measure the response of the flow to the symmetry of the forcing and χ = ∂S/∂θ will measure the corresponding susceptibility. Identification of relevant S can be made by studying the flow topology in the laminar and steady regime, at low Reynolds number, i.e. for Re 200.

The von Kármán flow topology: the laminar case

When θ = 0, the produced laminar flow is steady and composed of two toric recirculation cells separated by an azimuthal shear layer located at z = 0 (cf. figure 2(a)) reflecting the R π -symmetry of the system. Actually, each impeller pumps the fluid located near the rotation axis into its center before expelling it radially. This centrifugal pumping creates a toric recirculation cell in front of each impeller. In parallel, in each of the two toric cells, the fluid is rotating in the azimuthal direction following the nearest impeller. Then, when the motors frequencies differ, i.e. when θ = 0, the shear layer moves towards the slowest impeller, breaking the symmetry of the flow, but keeping the two-cells topology (cf. figures 2(b-d)). As θ tends to ±1, the flow continuously tends to a one-cell topology i.e. a topology with a single toric recirculation cell.

Identification of two order parameters and associated susceptibilities

We will use two different order parameters to quantify the distance of the flow to the R π -symmetry. First, we consider a local quantity: the z-position of the shear layer of the flow, i.e. the position z s of the stagnation point on the rotation axis ‡. The stagnation point position z s is only defined when the flow has the particular two-recirculation-cells axisymmetric structure and will not be suitable -because non-defined-to monitor the distance of the instantaneous flow to the R π -symmetry when it is turbulent. For this purpose, we prefer a second order parameter, a global quantity, the normalized and space-averaged angular momentum I(Re, θ, t) :

I(t) = 1 V V rdrdϕdz ru ϕ (t) π R 2 (f 1 + f 2 ) ( 1
)
where V is the volume of the flow. Practically, I(t) is computed from SPIVdata restricted to a meridian plane only but, since azimuthal flow fluctuations are ‡ Practically, for an axisymmetric flow, z s is measured as the axial position, ψ(r = 0, z s ) = 0, of the zero iso-surface of the stream function ψ(r, z) which is defined through (u r , 0, u z ) = ∇ × (r -1 ψ e ϕ ) in cylindrical coordinates. strong, time-average over several impeller rotation periods -statistically equivalent to spatial azimuthal averaging-estimate correctly the 3D-value of I(t) (see details in reference [25]).

Examples of time-series of I(t) in turbulent regimes are provided below in figure 3 and reveal that strong spontaneous symmetry fluctuations can occur (cf. figure 3(a) for Re = 126 000, θ = 0). We assume that ergodicity holds, meaning that the instantaneous turbulent flow is exploring along time its energy landscape according to its statistical probability. In this framework, the time average value I of I(t) is equivalent to a statistical mechanics ensemble average providing the average is performed over a long enough duration -comparable to the averaging time needed to measure z s -in order to correctly sample the slowest time-scales. The quantity I shows the major advantage of being, contrary to z s , a robust observable since it is defined for each instantaneous velocity field whatever the structure of the flow and therefore even for turbulent flows. When observing the statistics of the fluctuations of I, we will sometimes encounter non-Gaussian and even multi-peaked distributions. Therefore, we will consider the time average I, the most probable value(s) I * and the fluctuations I(t) of the global angular momentum I to characterize the mean, most-probable and instantaneous symmetry of the flow.

Then, using order parameter I, we define two susceptibilities of the flow to symmetry breaking χ I and χI , based respectively on the mean and most-probable values of I, as:

χ I = ∂I ∂θ θ=0 and χI = ∂I * ∂θ θ=0 . (2)
Another susceptibility may also be defined using z s as:

χ z = dz s dθ θ=0 . (3)
Finally, since our observations (cf. figure 4) have shown that I is strictly proportional to the mean altitude z s of the shear layer, with I/z s = 0.23±0.01 whatever Re, we may simply use the susceptibilities χ = χ I = 0.23χ z and χ = χI . In the sequel, we investigate the influence of turbulence on these order parameters -z s and I(t)and the associated susceptibilities as Re, and therefore turbulence intensity, increases from 10 2 to 10 6 .

About the statistics of the time-series and high-frequency noise

The statistics of the long time series of I(t) obtained from SPIV data represent a crucial step in our study but they are limited by storage space and computing time. Our results are partly based on the mean values I and its equivalent z s computed on mean velocity fields and these lowest moments are well converged for all the data presented. Since statistical physics generally considers the most probable value of a distribution rather than its mean value and since some signals appear multivalued (cf. figure 3(a) and figure 8 below), we are also interested in measuring absolute and relative most-probable values. However, the PDFs of the all raw signals I(t) are almost Gaussian and symmetric. Lowpass-filtering applied on I(t) with a typical 0.15Hz-cutoff removes some high-frequency noise, due in particular to the small scales of turbulence, and leads to PDFs with relevant non-Gaussian shapes. The best example is the high-resolution time-series corresponding to the highest fluctuation-level measured at Re = 39 000, presented in figure 8. The most-probable values are computed accordingly, but because of statistical limitations, these value are less converged than mean values.

Formal analogy with ferromagnetic systems

In the non-fluctuating laminar case, when θ = 0, I is strictly equal to zero due to the symmetry of the flow. In contrast, as θ drifts away from 0, the value of the angular momentum I becomes more and more remote from zero as the asymmetry of the flow grows. In the turbulent regimes, the same behavior is observed but with higher susceptibilities to symmetry breaking, i.e. higher sensitivity to θ.

In such a framework, the Reynolds number -or a function of it [26]-is the equivalent of a statistical temperature and we are entitled to propose the following formal analogy between ferromagnetic systems and von Kármán systems:

• Order parameter: magnetization M ⇔ angular momentum I;

• Symmetry breaking parameter: external applied field h ⇔ relative driving asymmetry θ;

• Control parameter: temperature T ⇔ Reynolds number Re, or a function of it.

3. Laminar and turbulent flow: influence of θ and Re

The laminar susceptibility

The evolutions of the symmetry parameters z s and I as a function of θ in the laminar flow at Re = 120 ± 25 are provided in figure 5. In this case I(t) = I(t) and z s are defined for each instantaneous field. We see that z s evolves almost linearly with θ from z s = -0.7 R (position of impeller 1) for θ = -1 to z s = 0.7 R (position of impeller 2) for θ = 1, through z s = 0 for θ = 0 (cf. figure 5(a)). Variation of I with θ is similar.

In the central region, the linearity is excellent and we measure χ z = 1.03 ± 0.08 and χ I = 0.240 ± 0.005.

Flow topology : the turbulent case

Increasing the Reynolds number, one expects to reach fully developed turbulence around Re = 10 000 [11]. In this turbulent regime and at θ = 0, the R π -symmetry is broken for the instantaneous flow. However, as usually observed for classical turbulence, this symmetry is restored for the time-averaged flow (cf. figure 6(a)). Then, as in the case of the laminar flow, when θ is varied, we observe the breaking of the R π -symmetry of the turbulent mean flow, the structure of which matches quite well the laminar flow topologies, however with a higher sensitivity to θ (cf. figures 6(b-d)).

The turbulent susceptibility

The main properties of z s and I with respect to θ can be observed in figure 7 which combines the measurements obtained at an intermediate Reynolds number Re = 67 000 ± 2 000 (red circles) as well as at the highest reachable Reynolds number in our experiment Re = 890 000 ± 140 000 (blue triangles). The left-hand-side figures give a global view over a large span in θ: at this scale the two high-Reynolds-number series appear very close to each other and show much higher susceptibilities than the laminar flow series (green squares): turbulence enhances dramatically the sensitivity of the flow to symmetry breaking. The turbulent susceptibilities are at least 10 times larger than the laminar susceptibility since z s saturates -the shear layer disappears into one impeller-for |θ| larger than a finite and small value θ c = 0.10 ± 0.02, close to the value predicted in [27] for flat disks and 10 times smaller than in the laminar case where the shear layer disappears asymptotically as θ → ±1 (cf. figure 2 and 5(a)). These disappearances of the shear layer into impeller 1 or impeller 2 correspond to the transitions of the flow topology from two recirculation cells -at small |θ|-to a single recirculation cell -at large |θ|. In figure 7(e), we observe two maxima for σ I , the level of fluctuations of I, around ±θ c , which are characteristic of these transitions. Such transitions have already been studied for other counterrotating von Kármán flows in subcritical cases [10,25].

The effect of the Reynolds number on turbulent flow responses is revealed by the ×10-zooms in θ on the right-hand-side of figure 7. It is clear that, in a narrow region |θ| 0.01, we encounter two major differences between the Re = 67 000 and Re = 890 000 experiments:

• The susceptibilities are very different: χ I = 42 ± 1 at Re = 67 000 and χ I = 9 ± 1 at Re = 890 000 (cf. figure 7(b,d)), respectively 175 and 37 times the laminar value.

• The fluctuation level σ I (θ) presents a sharp and narrow peak at θ = 0 for Re = 67 000, but not for Re = 890 000 (cf. figure 7(f)).

This peak is typical of all data in the intermediate Reynolds number range 20 000 Re 200 000 where the susceptibility is clearly higher than around the highest Reynolds number reached Re ∼ 10 6 . Note that outside the active region of high susceptibility and high fluctuations |θ| 0.01, there is no more Re-dependency for I(θ): in the medium range 0.01 |θ| 0.05, not too close from the two-cells/one-cell transition at ±θ c , all curves collapse on tilted plateaus which -if extrapolated-cross the θ = 0 axis at I 0 ±0.05 whatever Re.

Fluctuation level near θ = 0

For all Reynolds numbers in the intermediate range, we observe a sharp and narrow peak for the fluctuation level σ I (θ) around θ = 0. This is illustrated by two time-series at Re = 127 000, in figure 3(a) and (b), for θ 0 and θ = 0.0084 which correspond respectively to the top and the bottom of the fluctuation peak in figure 7(f). The difference between the two signals is striking and reveals that the fluctuation peak is due to the presence non-Gaussian intermittencies. These intermittencies have been observed in the range |θ| 5 × 10 -3 exclusively. The I(t) time-series at Re = 890 000 and θ = 0, in figure 3(c), corresponds however to the minimum fluctuation level of figure 7(f): the corresponding PDF is narrow, symmetric and has a single central maximum.

Mean vs. most-probable values for highly intermittent regimes

In the above paragraphs and figures we have presented the evolution of mean values of I(t) with θ for different Re. Globally, most-probable values differ from mean values only when the susceptibility is high and are, in absolute value, greater than mean values. This difference is especially noticeable for Re around 40 000. The data set which corresponds to the maximum fluctuation level is presented in figure 8 where I 0. It shows strong temporal intermittency between states of different angular momentum values revealed by the multi-peak PDF.

The most-probable values of I(t) are plotted in figure 9, with respect to θ very close to θ = 0. We observe that for Re = 37 000 ± 2 000, the slope χ for most-probable values is approximately twice the slope χ for mean values at θ = 0. We also reported points obtained from three high-resolution time-series at Re = 39 000 (figure 8) and 43 000. For these data sets, we plot the first two most-probable values ±I * of I(t) corresponding to the two almost symmetric peaks of the PDFs with are exchanging their relative importance around θ = 0. These values are significantly higher than those of the Re = 37 000-data series. This can be due to either the proximity of the critical value Re c or to the limited statistics of the Re = 37 000-data series. Anyhow, these results suggest that the corresponding susceptibility χ is diverging at Re = 41 000 ± 2 000 as it will be discussed in the next section.

Divergence of the susceptibility at Re 40 000

From the above results it is clear that the effect of the Reynolds number is to be studied for θ = 0. The present section deals with the susceptibilities and the fluctuation level at θ = 0 and for Re covering the whole study range 10 2 -10 6 .

In figure 10(a), we plot susceptibility measurements -χ and χ, cf. equation 2-as a function of Re. Figure 10(b) shows the inverses of these values and figure 10(c) the corresponding standard deviation σ I of I(t). The susceptibility χ based on the most probable values I * of I(t) and the fluctuations of I(t) are both maximum for Re 40 000. However, as reported earlier [21], the susceptibility χ based on the mean values of I(t) reaches a maximum for higher Re 90 000. We also note that χ χ, whatever Re, and that the highest measured susceptibilities are of the order of the highest measurable values considering the θ-precision of our setup.

All the above results suggest a divergence of χ(Re) which clearly appears with the fits in figures 10(a) and (b): χ(Re) diverges at Re = Re c = 40 000 ± 5 000 with critical exponent -1 [illustrated in figure 10 with respect to 1/ log(Re) (see discussion in section 6.2)] and resembles a classical magnetic susceptibility divergence at the ferromagnetic/paramagnetic phase transition. In the next section, we investigate further this transition by looking at the equivalent of the spontaneous magnetization of magnets, i.e. how, around Re c , the angular momentum I(t) spontaneously and dynamically transits between different finite values. Close to the susceptibility divergence, we observe a complex dynamics for the global angular momentum I(t), as illustrated in figure 8 at Re = 39 000 and θ = 0. Indeed, one observes that I(t) does not just fluctuate randomly around zero -its mean valuebut shows a tendency to lock for some times, preferentially around ±I * with I * of the order of 0.03-0.04, estimated from the main peaks of the PDF of the filtered signal I f (t) (cf. figure 8(b)).

In order to understand the nature of the fluctuations of I(t), we have first checked that I(t) is uncorrelated with the control parameter θ(t) (correlation coefficient C(I, θ) = I(t) θ(t)/σ I σ θ < 0.02, cf. figures 8(a,c)) and thus reveals only spontaneous fluctuations. Actually, I(t) fluctuates with two separate time scales: (i) fast fluctuations related to "traditional" small-scale turbulence and (ii) intermittencies corresponding to residence time of typically few tens of seconds and certainly due to the proximity of the susceptibility divergence. If one performs a time-average of the velocity field over one of these intermittent periods, one obtains a time localized "mean" flow, analogous to what is obtained for true mean flows when θ = 0 as presented in figures 6(b-c), and with its own level of spontaneous symmetry breaking. To go one step further with the correspondence between I(t) and the amount of symmetry breaking or the position of the shear layer, we have computed bin-averages of the velocity field for a time-series Figure 11. Bin-averaging of the velocity field with respect to I(t) for Re = 39 000 and θ = 0 (same data series as in figure 8). Three bin-averaged velocity fields are presented for bins: (a) -0.075 < I f (t) < -0.050, (b) -0.0125 < I f (t) < 0.0125 and (c) 0.0375 < I f (t) < 0.0625. (d) Corresponding stagnation point axial position z s measured for 11 such bin-averaged fields for I between -0.075 and +0.075 (green stars), superimposed to regular time-averaged data for varying θ at Re = 67 000 (red circles) and Re = 890 000 (blue triangles) as in figure 4(b). I(t) in 11 bins ranging from -0.075 to +0.075. In figures 11(a-c), we have plotted three examples of such bin-averaged fields. The pattern is R π -symmetric for the bin around zero (cf. figure 11(b)) and non-symmetric otherwise (cf. figures 11(a) and (c)). From the resulting series of conditionally-averaged velocity fields we measure the stagnation point axial position z s and report it as a function of the mean bin value (cf. figure 11(d)). There is a very good matching between these bin-averaged data and the timeaveraged data already presented in figure 4 for different θ: beyond the proportionality between z s and I, I(t) can thus definitively be used as a quantitative instantaneous symmetry measurement. For example, the two symmetrical most-probable values ±I * of I(t) correspond typically to |z s | 0.2 R.

From the above observations, we can conclude that, very close to Re c , the turbulent flow explores a band of metastable symmetry breaking patterns evidenced by -I * I f (t) +I * , where I f (t) is the 0.15Hz low-pass filtered value of I(t), the three most-visited states being -I * , 0, and +I * . Farther from Re c , signals and PDFs also reveal coexisting preferred values for I(t) when θ is very small.

The "momentization"

The dynamical spontaneous symmetry breaking observed in our experiment can be seen as the coexistence of multistable/metastable steady states of different symmetry which are explored dynamically along time. By analogy with magnetization we propose to describe this as successive "momentization" patterns observed in time, thanks to ergodicity which arises from the very high level of fluctuations. With respect to magnetism, the main difference is that the spontaneous momentization occurs along time and mainly very close to the critical point. Furthermore, owing to the statistical limitations, we cannot conclude if momentization occurs only on one side of the transition -as it is for the ferro/para transition-or on both sides.

Statistics of the dynamic multistability

The metastable/multistable states of finite lifetime (τ ) observed in time-series such as in figure 8(a) can be characterized by their angular momentum averaged over their lifetime and noted I τ , to avoid confusion with the full-signal mean value I. The I τ are distributed along the angular momentum axis and we note I(Re, θ) the I-interval containing the distribution for experiments at (Re, θ). Multistability is systematically observed in the high susceptibility zone limited by Re ∈ I Re ∼ [20 000, 200 000] and θ ∈ I θ ∼ [-0.01, +0.01]. Whatever Re and θ in I Re and I θ , we observe that I(Re, θ) is confined in [-I 0 , +I 0 ] (I 0 has been defined at the end of section 3.3 and is shown in figure 7). When Re Re c , the interval I(Re, θ) [-I * , +I * ] is symmetric and close to [-I 0 , +I 0 ] (cf. figure 8): the whole band is almost populated. However, farther from Re c , the interval I(Re, θ) is smaller and evolves inside [-I 0 , +I 0 ] with θ.

Complementary investigation about the azimuthal structure of the flow

Before discussing our results in the next section, it is important to note that all the above results concern 2D-velocity-fields in a meridian plane: the azimuthal spatial structure is ignored by our 2D-SPIV process since the azimuthal patterns -moving and fluctuating-are averaged along time and not resolved in space at all. In this section, we investigate some preliminary experimental results about the 3D-structure of the flow, gathered by visual observation with bubble seeding and white light [11,25]. With a high density of bubbles we mostly see the patterns on the outer cylinder which are the radial vortices due to the shear-layer destabilization. In practice, we experience that, even in the flow is turbulent, the azimuthal number m of vortices can be a well defined integer over some ranges of Reynolds number and/or some periods of time. So, exactly as for small dynamical systems, changes in m can reveal an Eckhaus instability for the 3D-mean pattern §.

This dynamics has been observed for 22 000 Re 300 000. Despite the vortices are strongly fluctuating, we can report the following observations of the average number m of vortices around the perimeter of the shear layer at θ = 0:

• m = 3 for 22 000 Re 47 000,

• m = 4 for 86 000 Re 300 000.

Between these two regions, the fluctuations are generally too strong to draw precise conclusions, but at least a clear observation of m = 4 has been made once at Re 51 000. So, we conclude that some kind of m = 3 ↔ m = 4 Eckhaus transition occurs between Re = 47 000 ∼ Re c and Re = 86 000 ∼ 2Re c , in the region of very high susceptibility.

Discussion

We have reported the experimental study, for a von Kármán swirling flow, of (i) the response of the mean flow to a continuous breaking of the R π -symmetry of the system and (ii) the spontaneous symmetry fluctuations of the instantaneous flow, from the laminar regime at Re ∼ 10 2 to the highly turbulent regime at Re ∼ 10 6 . The divergence of the susceptibility χ(Re) for Re = Re c = 40 000 ± 5 000 reveals the existence of a phase transition or a bifurcation. In the following, we discuss the implications of this phase transition for turbulence, the critical behaviors and the analogy with the paraferromagnetic transition.

A Reynolds-dependent turbulence

This new transition makes some non-dimensional characteristic quantities of the flow still Re-dependent at very high Re. The turbulence is thus definitively not fully developed in this closed turbulent flow, contrary to the most common observations in this Reynolds number range. In other closed flows, the literature reports some examples of transitions occurring at Reynolds number of the order of 10 5 , a decade in Re beyond where the turbulence already looks like fully developed. First, in various type of turbulent fluid experiments, several authors report transitions and symmetry breaking at high Reynolds number (cf., e.g., references [START_REF] Chillá | [END_REF]5,6,7]). Moreover, in a similar von Kármán flow, de la Torre and Burguete [12,13] observed a bistability between two broken O(2)-symmetry states with z s /R = ±0.19. Although the slow dynamics is very different from what is reported here, both experiment are likely to belong to a common framework (see below). § We use here the word "mean" not for the 2D-or axisymmetric-time-average, as everywhere in the present article: here, "3D-mean" involves time averaging in a rotating frame following the azimuthal pattern. Such image processing will be a challenge, implying movies of the whole system [28] or full 3D velocimetry.

Finally, in a liquid helium von Kármán flow, Tabeling et al. report a local peak in the flatness of velocity derivative around R λ = 700 -corresponding to Re = 2 × 10 5 [8]. An interpretation in terms of a second-order phase transition was proposed [9]. Velocity derivatives cannot be measured in our experiment because of the low spatial resolution of our SPIV system. However, if one could establish that the observations of Tabeling et al. are indeed connected with ours, then our work would prove that their transition is related to the global structure of the mean flow and thus does not result from the breakdown of small-scale vortical structures as proposed in reference 6.2. Critical behaviors and analogy with magnetism 6.2.1. The control parameter: a function of log Re -Since our study covers 4 orders of magnitude in Reynolds number, it is natural to acquire and plot results along a logarithmic scale. However, we could search for critical behaviors with either Re or log Re as control parameter. Indeed, critical behaviors are asymptotic behaviors around a given threshold Re c and therefore both (Re -Re c ) and (log Re -log Re c ) are mathematically equivalent. However, it is very common in non-linear physics, that asymptotic behaviors or normal forms correctly describe the dynamics, not only asymptotically, but over finite -and even large-ranges around thresholds using one particular control parameter. This is exactly what happens here and because of the quality of the 1/ χ(Re) fit in figure 10(b), we conclude that a relevant control parameter for this critical phenomenon can be log Re or a gentle function of it.

Within the formal turbulence/magnetism analogy, we actually use Castaing's proposition [26] for a temperature of turbulent flows which is T ∼ 1/ log Re.

Critical behavior of the susceptibility -

The divergence appears clearly in figure 10(b) where 1/ χ(Re) goes to zero at Re = Re c = 40 000 ± 5 000 in a way that can be fitted with a -1 critical exponent by: χ

∝ 1 log Re - 1 log Re c -1 (4)
with different prefactors below and above Re c . This expression can be interpreted as the combination of the Castaing statistical temperature of turbulence T ∼ 1/ log Re with the classical mean-field-theory expression for susceptibility at a critical point:

χ ∝ |T -T c | -1 ,
Expression (4) has already been successfully used in reference [21] to describe the "mean" susceptibility χ computed with mean values -instead of most probable values for χ-of I(t) and with a more limited set of data. As already underlined, this rather fair agreement between mean-field-theory and our turbulent measurements is at first surprising. However, recent theoretical developments suggest that in turbulent von Kármán flow the interactions are long-range like [29,30] supporting therefore the idea that mean field theory is nearly applicable in our system.

Momentization vs. magnetization

Within the proposed analogy between our observations and the para-ferromagnetism transition, the equivalence between the magnetization and the momentization is the most questioning aspect. In particular, the shape of the I(θ)-curves is more similar to the magnetization response curves of a single magnetic domains than to poly-domains.

For the para-ferromagnetism transition, the magnetization increases with the distance below the critical temperature while here, as stated in section 4.3, the spontaneous momentization is finite, limited in amplitude by I 0 , and decreases with the distance to Re c . The system may share some characteristics with reentrant noise-induced phase transition similar to that observed in the annealed Ising model [31,32]. Our system is definitively of finite size and it is thus coherent that both the level of fluctuation σ I (Re c ) and the band of spontaneous fluctuations for I(t) -which corresponds to a finite excursion of the shear layer position along z-remain finite at critical point (figure 10(c)).

In their von Kármán experiment at θ = 0 and Re 4 × 10 5 , de la Torre and Burguete ([12, 13]) also observe multistability of three and only three states with can be interpreted, within our framework, as I τ = -I * , 0 and +I * . While the amplitude of the I τ = 0 mode is constant along time, their system switches between the two I τ = ±I * modes with a very slow dynamics (hours) and very fast transitions. This dynamics does not seem to depend on Re once Re 2 × 10 4 , the symmetric I τ = 0 state being dominant below. Within our framework, we could propose that, in this experiment, the susceptibility is infinite whatever Re above 2 × 10 4 (cf. figure 5(a) of reference [12]).

Specificity of turbulent fluctuations with respect to classical thermal noise

In our turbulent system, the observed time-dynamics seems to be due to the very high level of intrinsic turbulent fluctuations, which are very different from classical thermal fluctuations that pilot the dynamics of magnetic momentums. These turbulent fluctuations are considered as the noise or the statistical temperature origin of our system. In classical close-to-equilibrium thermodynamical systems, thermal fluctuation are generally considered as regular additive noise exactly as for the classical treatment of instabilities in the presence of noise where an external noise -additive or multiplicative-is introduced as a perturbation of a regular system [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF]. For turbulent flows in general and for our von Kármán flow in particular, the situation is different: the flow itself is intrinsically highly fluctuating so that fluctuations are of the same order of magnitude than the mean values [25] and cannot be treated as perturbations. So, the specific problem our experiment addresses is the behavior of the instability of an average global quantity such as I(t) in an intrinsically fluctuating or noisy system. This is similar to the problem of the experimental dynamo instability [START_REF] Monchaux | [END_REF]35] which can only occur in turbulent flows because of the nature of the existing conducting fluids of the Universe. Therefore, non-classical behaviors are likely to be observed for bifurcations or phase-transitions in such highly-turbulent flows and our system represents a unique tool to study these transitions.

Figure 1 .

 1 Figure 1. Schematic view of (a) the experimental setup and (b) the impellers blade profile. The arrow indicates the rotation sense. (c) 3D-view and symmetry : the system is symmetric with respect to any R π -rotation of angle π around any line in the equatorial plane which crosses the rotation axis.

Figure 2 .

 2 Figure 2. Maps of velocity fields in the meridian plane (r, z) -the rotation axis is vertical at r = 0 and the r ↔ -r symmetry of the maps reveals the flow axisymmetryfor the laminar (Re = 120 ± 25) von Kármán flow for different values of the symmetry control parameter θ: (a) θ = 0, (b) θ = -0.11, (c) θ = -0.19, (d) θ = -1. The color maps, from blue to red ("jet" colormap), the azimuthal velocity u ϕ , whereas the arrows map the (u r , u z ) field. The resolution of the fields has been reduced by a factor 2 for better visibility. Positions r and z are given in units of vessel radius R.

Figure 3 .

 3 Figure 3. Global angular momentum I(t) as a function of non-dimensional time f.t for experiments performed at Re = 126 700 ± 200 with (a) θ = 0 ± 10 -4 and (b) θ = -0.0084 and (c) at Re = 890 000 and θ = 0. Green thin signals are SPIV data I(t) sampled at 15 Hz and blue thick signals correspond to low-pass filtered data I f (t). Corresponding right figures shows the probability density functions (PDF) of the signals with same color code. The red lines indicate the most-probable values of I f (t).

Figure 4 .

 4 Figure 4. Experimental observation of proportionality between time-averaged global angular momentum I and mean shear layer altitude z s for (a) laminar flow Re = 120 ± 25 (green squares, I/z s = 0.235 ± 0.005) and (b) turbulent flows at Re = 67 000 (red circles, I/z s = 0.227±0.005) and Re = 890 000 (blue triangles, I/z s = 0.23±0.01). Each data series is obtained through plotting I(θ, Re) as a function of z s (θ, Re) for varying θ and constant Re. At high Re, this proportionality holds for -0.1 θ 0.1, i.e. while the shear layer position z s is not saturated to the impellers positions ±0.7 R.

Figure 5 .

 5 Figure 5. (a) axial position z s of the stagnation point and (b) global angular momentum I = I as a function of θ for the laminar flows at Re = 120 ± 25.

Figure 6 .

 6 Figure 6. Maps of mean velocity fields of the turbulent von Kármán flow at Re = 890 000 for different values of θ: (a) θ = 0, (b) θ = -0.0036, (c) θ = -0.0147, (d) θ = -1, with same layout as figure 2. The r ↔ -r symmetry of the maps reveals that the time-averaged mean fields are axisymmetric.

Figure 7 .

 7 Figure 7. Turbulent regimes at Re = 67 000 ± 2 000 (red circles) and Re = 890 000 ± 140 000 (blue triangles): axial position of the stagnation point z s (a-b) and global angular momentum I -mean I(t) (c-d) and standard error σ I (e-f)-as a function of θ. Wide span -0.3 < θ < 0.3 on the left side (a,c,e) and ×10 zoom -0.03 < θ < 0.03 on the right side (b,d,f). Data for the laminar regime at Re = 120 (green squares) from figure 5 are also plotted for comparison. The data have been symmetrized with respect to θ = 0. Dashed lines in (d) indicates the tilted plateaus independent of Re and their extrapolations to ±I 0 .

Figure 8 .

 8 Figure 8. (a) Global angular momentum I(t) for an experiment performed at Re = 39 000 for θ = θ(t) = 0. The green thin line is SPIV data I(t) sampled at 15 Hz and the blue thick line corresponds to 0.15 Hz low-pass filtered data I f (t). (b) Probability density functions (PDF) with the same color code. Blue dashed lines indicate the two most-probable values ±I * of I f (t) and green dash-dot lines indicate the mean value I. (c) Time-series of θ(t) as produced by the motor regulation in unit of 10 -3 : the correlation with I(t) is negligible (see text).

Figure 9 .

 9 Figure 9. Mean I (red circles) and most-probable I * (blue stars) values of the global angular momentum I(t) for Re ∼ 40 000. Close symbols come from a set of moderateresolution data at Re = 37 000±2 000 (f = 12Hz, 1 200 SPIV fields acquired at 1.7Hz), whereas open symbols come from high-resolution time-series at Re = 41 000 ± 2 000 (f = 8Hz, 9 000 or 27 000 SPIV fields acquired at 15Hz as, e.g., in figure 8) acquired in a different glycerol-water mixture. The data have been with respect to θ = 0 and, for the data point at θ 0 of the second series, we plot the two most-probable values ∼ ±I * corresponding to the two peaks of the PDFs as in figure 8(b).

Figure 10 .

 10 Figure 10. Reynolds number dependence of the symmetry-breaking susceptibility of the von Kármán flow at θ = 0 (a) and its inverse (b). Black open symbols correspond to susceptibility χ = χ I = 0.23χ z obtained from mean values I (squares) and z s (diamond); red full circles stands for χI obtained from most-probable values I * of I. The red full star corresponds to the region (Re = 41000 ± 2000) where infinite susceptibility is inferred from the strong multistability of the signal (cf. figures 8 and 9). Dashed lines are fits of the most-probables values above and below the star as linear functions of 1/ log(Re) as proposed in the discussion of section 6.2. (c) Standard deviation σ I of the unfiltered global angular momentum I(t) at θ = 0.

Acknowledgments

We thank Séverine Atis and Lise Divaret for their participation to the data acquisition. We thank Cécile Wiertel-Gasquet for the interfacing and control of the experiment. We thank Kirone Mallick, Sébastien Aumaître, Javier Burguete and Frédéric Moisy for fruitful discussions. PPC was supported by Triangle de la Physique and EH by ANR SHREK (ANR-09-BLAN-0094-03).


P.-P Cortet

D Bonamy

F Daviaud

O Dauchot

B Dubrulle

M Renouf

Relevance of visco-plastic theory in a multi-directional inhomogeneous granular flow

Keywords: 45, 70, -n -Granular systems 83, 80, Fg -Granular solids 45, 05, +x -General theory of classical mechanics of discrete systems

We confront a recent visco-plastic description of dense granular flows (Jop P. et al., Nature, 441 (2006) 727) with multi-directional inhomogeneous steady flows observed in nonsmooth contact dynamics simulations of 2D half-filled rotating drums. Special attention is paid to check separately the two underlying fundamental statements into which the considered theory can be recast, namely i) a single relation between the invariants of stress and strain rate tensors and ii) the alignment between these tensors. Interestingly, the first prediction is fairly well verified over more than four decades of small strain rate, from the surface rapid flow to the quasi-static creep phase, where it is usually believed to fail because of jamming. On the other hand, the alignment between stress and strain rate tensors is shown to fail over the whole flow, what yields an apparent violation of the visco-plastic rheology when applied without care. In the quasi-static phase, the particularly large misalignment is conjectured to be related to transient dilatancy effects.

Despite granular flows ubiquity in nature and industry, their description still remains an open issue. Three regimes are usually distinguished: The rapid dilute flow -gaseous-like-regime, where grains interact through binary collisions, is described within the framework of the kinetic theory [1]; the slow flow -solid-like-regime, where only friction between grains is relevant, is described using soil mechanics [START_REF] Nedderman | Statics and Kinematics of Granular Materials[END_REF]. In between, there is a denseflow -liquid-like-regime, where inertia becomes important, friction remaining however relevant. This last regime has been widely investigated experimentally and numerically (see [START_REF] Midi | [END_REF] for a review) in various flow configurations and has been the focus of many theoretical works [4][5][6][7][8]. However, while these descriptions reproduce some of the features observed experimentally, a unified description is still missing. In particular, the possibility for a single and local constitutive law within the framework of continuum mechanics is still debated.

Progresses have been made recently in this context. Dimensional analysis [START_REF] Midi | [END_REF]9,10] was used to show that, in simple incompressible unidirectional uniformly sheared flows, bead diameter d, grain density ρ, global shear stress τ , global pressure P and global shear rate γ should be related as τ /P = µ(I), where µ(I) is an effective friction coefficient, function of the so-called Inertial Number I = γd/ P/ρ expected to control locally the rheology of the material. Such relation has first been evidenced, numerically, in plane shear [9,10] and, experimentally, in inclined plane [START_REF] Midi | [END_REF] configurations. Local tensorial extension of this relation was recently proposed by Jop et al. [11] as a constitutive law for dense granular flows:

τ ij = µ(I) P | γd | γd ij with I = | γd |d P/ρ , (1)
where

τ ij = σ ij + P δ ij denotes the deviatoric part of the stress tensor, γij = 1/2 (∂ i v j + ∂ j v i ) is the strain rate tensor, | γd | = 1 2 γd ij γd ij the norm of its deviatoric part γd ij = γij -1/2 γkk δ ij , P = -1/2
σ kk the local pressure and v i the components of the velocity. In the compressible case, a consistent local rheology would imply in addition a univocal dependance ν = ν(I) of the volume fraction on the inertial number [10]. Introducing eq. ( 1) in the momentum equation, Jop et al. succeeded to fit the surface velocity profile for a steady unidirectional flow down an inclined plane with walls. More recently, relation (1) was confronted with simulations of complex multi-directional non-steady inhomogeneous flow -the collapse of a cylinder of granular matter onto a horizontal plane [12]-and was shown to hold fairly well within the fast flowing region i.e. for I ranging from 10 -2 to 0.8. In this last work, relation (1) has however been verified only for positions where stress and strain rate tensors had parallel principal directions, excluding 10% of data from the test. Therefore, despite these various experimental and numerical results, the relevance and limits of the visco-plastic relation (1) remain to be investigated in granular flows involving small values of inertial number I -in which the stress and strain rate tensors are particularly prone to be non-aligned-as e.g. avalanches flows where large scale non-local effects are expected [13][14][15][16].

In the present letter, we confront this visco-plastic description of granular flows with steady surface flows observed in a 2D rotating drum, using extensive numerical simulations providing large statistics. The motivation for studying this specific flow is double since, still remaining in a "simple" 2D geometry, i) it provides a general highly multi-directional inhomogeneous situation and ii) it offers a wide range of inertial number: 10 -6 < I < 5 × 10 -2 . In such a flow, we report an apparent failure of the constitutive relation (1) over the whole range of I, when applied in the "natural" frame of the flow. To overcome this failure, we highlight the fact that eq. ( 1) is actually composed of two independent results, namely a) a single relation between the invariants of stress and strain rate tensors and b) the alignment between these tensors, that constitute its two underlying fundamental hypotheses. The system of these two conditions actually recasts eq. ( 1) in a more fundamental formulation. Then, the norms of the stress and strain rate tensors are found to be univocally related through the invariant relation derived from eq. ( 1) over the whole range of I, meaning condition a) is verified over the whole flow. Finally, the failure of the tensorial relation ( 1) is shown to result from a misalignment between the stress and strain rate tensors -meaning condition b) fails-contrary to what is assumed in previous works. Possible origins of this misalignment are discussed: in the deep "solid" part of the flow, the total deformation experienced by the granular packing during its lifetime is less than one-bead diameter. This induces transient dilatancy effects -from a Lagrangian point of view-that could be at the origin of the persistence of the observed large misalignment.

The simulations have been performed using the Non-Smooth Contact Dynamics method [17]. The scheme has been described in detail in [18] and is briefly recalled below. An immobile drum of diameter D 0 = 45 cm is halffilled with 7183 rigid disks of density ρ = 2.7 g cm -3 and diameter uniformly distributed between 3 and 3.6 mm. The weak polydispersity introduced in the packing prevents 2D ordering effects. The normal restitution coefficient between two disks (respectively, between disks and drum) is set to 0.46 (respectively, 0.46) and the friction coefficient to 0.4 (respectively, 0.95). Once the initial packing is stabilized, a constant rotation velocity Ω ranging from 2 to 12 rpm is imposed to the drum. After one round, a steady continuous surface flow has developed. Then, one starts to capture 400 snapshots equally distributed over one rotation of the drum. For each bead, one records the position r of its center of mass and its instantaneous velocity c over the time step δt = 6 × 10 -3 s. For each rotation velocity, we performed 20 experiments starting from different initial packings. The reference frame is defined as the frame rotating at Ω with the drum that coïncides with the reference frame 0 = ( e x , e z ) (cf. fig. 1), fixed in the laboratory, where e x (respectively, e z ) is parallel (respectively, perpendicular) to the mean flow free surface. The mean values of the various Eulerian fields (volume fraction, velocity and stress tensor) are obtained, at the grain scale (see [18]) and without any further coarse-graining, after both time and ensemble averaging over all the 400 frames of all the 20 experiments.

Steady surface flows in rotating drums provide a continuum of flow regimes from "quasi-static" to "liquid". Figures 1(b) and (c) present typical velocity profiles v x (x = 0, z) measured in our numerical simulations "at the center" of the drum (i.e. along the z -aligned slice at x = 0). Two phases can be identified: A flowing layer exhibiting a linear velocity profile and a "solid" phase experiencing creep motion where v x decays exponentially with z. Such velocity profile corresponds well to experimental measurements in rotating drum [START_REF] Midi | [END_REF][19][20][21] as well as in heap geometry [START_REF] Midi | [END_REF]22]. It is worth mentioning that the characteristic length ξ associated with the exponential creep is found to be ξ = 12 ± 2 in our numerical drum over the range of explored rotating velocity. This is significantly larger than that reported by Crassous et al. [23] in heap geometry, ξ = 1 ± 0.2. This difference may emphasize the role played by the forcing on the selection of the velocity profiles within the bulk of the material. Indeed, the forcing in rotating drum and heap are very different since, from a Lagrangian viewpoint, the gravity is rotating in the first case while it remains frozen in the second. The volume fraction is found to be almost constant, around 0.827 ± 0.05, in the creeping layer while it decreases slightly in the flowing layer, down to 0.80 at z = -4 (cf. fig. 1(d)) as expected from dilatancy effects.

The full 2D fields of deviatoric stress τ ij (x, z), deviatoric strain rate γd ij (x, z) and pressure P (x, z) are then computed to test the relevance of constitutive equation (1). Figures 2(a 1) even does not work well (though better) when restricting data to the slice x = 0 located at the center of the drum, where the flow is expected to be locally unidirectional (// to e x ) and almost homogeneous (∂ x 0) as in the plane shear configuration. Moreover, the shape of these curves does not correspond to the logarithmic shape previously evidenced in [START_REF] Midi | [END_REF]10,11], especially, because of the unexpected bump of τ xz | γd |/ γ d xz P for I ranging between 10 -4 and 10 -5 . The multi-directionality of the considered flow implies that no reference frame appears more natural than another to study tensors. Therefore, it is reasonable to recast tensorial equations into reference frame invariant relations. For granular flows, searching for a rheology consists mainly in relating the deviatoric stress tensor τ ij to the deviatoric strain rate tensor γd ij . Since both of these tensors are symmetric with zero trace, both of them are fully characterized through one invariant only, e.g. their norm, and their principal directions. Accordingly, the rheology of eq. ( 1) can be recast in the following invariant way: where |τ | = 1 2 τ ij τ ij is the deviatoric stress tensor norm. In unidirectional flows, such as those observed in plane shear geometry [10] or down to rough inclines [11], condition b) is naturally ensured by the flow symmetries so that verifying eq. ( 1) reduces to verifying point a) only. As we will see thereafter, this is no more the case for a multidirectional flow for which the alignment of the stress and strain rate tensors cannot be assumed [24,25]. The split of eq. ( 1) into the two independent conditions a) and b) recasts the rheology of Jop et al. into its two underlying fundamental assumptions which are completely independent. In the sequel, we will see that the invariant form of the visco-plastic rheology -condition a)-actually holds well over the whole rotating drum and the whole range of I either condition b) holds or not. This reformulation of eq. ( 1) actually allows to test the visco-plastic theory for very small values of I down to 10 -5 where the condition on stress and strain rate tensor alignment is largely violated.

Let us first consider point a). In fig. 3, we have plotted the invariant effective friction coefficient |τ |/P as a function of the inertial number for the same data as in fig. 2. We see that all the data appear to collapse quite well on a single master curve µ(I) independently of the rotation velocity Ω. The shape of this master curve and the value of the corresponding friction coefficient are perfectly consistent with the ones already evidenced (µ = 0.30 ± 0.10 for I < 0.1) in different granular unidirectional flows [START_REF] Midi | [END_REF]10,11]. It is important to emphasize that this collapse is a non trivial result since we have plotted data for the whole 2D highly inhomogeneous fields. In figs. 3(c) and (d), we restrict data to the slice x = 0 located at the center of the drum where the flow is expected to be almost "plane shear like". We see that the |τ |/P profiles for the three considered rotation velocities fall with a very good precision on the same master curve µ(I, x = 0) for values of I as small as 10 -5 . Despite this last remark, this collapse is a strong evidence for the possible relevance of the constitutive law of eq. ( 1) in steady inhomogeneous multi-directional flows. More importantly, it is the first time such frictional rheology seems to hold over such a low inertial number range (10 -5 < I < 5 × 10 -2 ) whereas it is often expected to fail because of the jamming transition for I < 10 -3 .

We already said that in compressive flows a local rheology would imply a single dependence of the volume fraction ν with the inertial number I. In fig. 4(a), we see that the data corresponding to the whole 2D flow collapse indeed quite well on a master curve ν = ν(I) supporting partially the proposed rheological approach.

The second step is now to check point b) i.e. the collinearity between the deviatoric stress τ ij and strain rate γd ij tensors. Figure 4(b) shows the angle difference δθ between the principal directions of τ ij and γd ij as a function of the inertial number. The collinearity is not verified since δθ ranges from ±10 • , for the larger inertial numbers (5 × 10 -4 < I < 5 × 10 -2 ) corresponding to the "liquid" surface flow, up to 35 • in the deeper part of the flow, corresponding to smaller inertial numbers (10 -6 < I < 5 × 10 -4 ). Such increase in the dispersion of the tensors misalignment, as I decreases, has also been reported in shear zones [24]. More importantly, there is no univocal relation ruled by the inertial number such as δθ = δθ(I). Consequently, the fact the tensorial law ( 1) is not verified contrary to its invariant version, relation a), on the norms of the tensors holds is the direct consequence of the misalignment of the stress and strain rate tensors. Moreover, it shows, looking at the larger inertial numbers in figs. 2 and 4(b), that a 10 • discrepancy between principal directions is already significant to make the visco-plastic tensorial rheology, in its non-invariant version of eq. ( 1), flatter.

In fig. 3, the collapse of the data on a master curve µ(I) is not perfect and one may wonder if the dispersion around the "mean" friction coefficient µ(I, x = 0) measured at the center of the drum is only due to a lack of statistics. Therefore, in fig. 5(a), we plot the map of the deviation of the effective friction coefficient from the one at the center of the drum. One can observe a spatial structuration with a quadripolar pattern. This spatial dependence of the profile |τ |/P = µ(I, x) with the position x may be the consequence of non-local effects reflecting the long range influence of the boundaries through force-chains, but remains however an unexplained feature. It is worth noticing that a pure geometrical dependence of |τ |/P has been reported in [25] where the interest of studying the present regime of mixed geometrical and shear rate dependence was highlighted. In parallel, in figs. 5(b) and 6(a), we see that the cloud of data δθ vs. I shows as well a non-trivial spatial structure. There is however no evidence for any link between the spatial structures of the angle δθ and of the deviation around µ(I, x = 0), the origins of which are still to be understood.

Previous studies [24,25] have shown that the stress and the strain rate tensor principal directions can differ significantly due to transient dilatancy effects as long as the initial conditions of the packing have not been forgotten, i.e. as long as the local deformation is less than ∼ 1. Even if our flow is stationary, the transient effect reported above could be the explanation for the large misalignment observed in the deep part of the rotating drum (cf. fig. 6(a)). Indeed, in the "solid" part of the granular matter, a given "Lagrangian" packing is created at the bottom of the surface avalanche. Then, it has a short lifetime reduced to half a drum rotation period 1/2Ω corresponding to the time needed for him to reach the top of the avalanche where it is destructed. During this process, the deformation experienced by the packing is quite small: in fig. 6(b), we plot the map of the logarithm of the norm of the strain rate tensor times the half drum rotation duration log | γd |/2Ω as well as a black line delineating where the deformation δ = | γd |/2Ω equals unity. The "solid" granular matter under this line is experiencing a transient deformation since it is less than one. Finally, we see that this "transient" zone matches quite well the zone (cf. fig. 6(a)) where |δθ| is larger than 10 • making transient dilatancy effects a relevant explanation for the persistence of large alignment discrepancies. This observation does not imply that the spatial structure observed for the angle δθ is due to transient effects. On the contrary, this structure is created by the specific rotating-drum boundary conditions, the small deformations experienced by the matter explaining then the fact that this structure of δθ(x, z) = 0 can live on.

Conclusion. -In this letter, we have confronted a recent visco-plastic theory for dense granular flow to numerical steady flows in rotating drums. The existence of a univocal relation between a scalar local friction coefficient -defined as the ratio between the norm of the deviatoric stress tensor to the pressure-and the socalled inertial number conjectured in the theory is in first approximation fairly well verified over the whole drum, in both the creeping "solid" phase and the avalanching phase at the free surface, for inertial numbers ranging from 10 -5 to 5 × 10 -2 . Small dilatancies evidenced at the flow surface are also found to be univocally related to the inertial number. On the other hand, the local alignment between deviatoric stress and strain rate tensors prescribed in the visco-plastic theory is found to fail significantly over the whole range of inertial number. This leads us to argue that the visco-plastic rheology proposed in [11] is extremely efficient to describe quasiunidirectional flows as those investigated in [9][10][11] but remains unapplicable, when applied in its non-invariant form (eq. ( 1)), to highly multi-directional flows such as those observed in rotating drum. On the other hand, it is particularly remarkable that the Jop et al. rheology, when considered in its invariant form, holds for very small values of I down to 10 -5 and that even if the stress and strain rate tensor misalignment is large. This result suggests in particular to reconsider the previous apparent failure of the visco-plastic rheology observed for small I in other flow geometries [10,25,26] which could have been a partially wrong conclusion. The origin of the observed misalignment between tensors is conjectured to be the consequence of transient dilatancy effects linked to the small deformation experienced by the "solid" part of the granular matter when considered from a Lagrangian point of view. Finally, a quadripolar spatial structuration of the dispersion around the perfect collapse of the data on the visco-plastic invariant rheology is tentatively attributed to non-local effects set by boundary conditions. In this context, stress distribution in immobile granular packings was found to be well described by Linear Elasticity [27]. It will then be interesting to see whether the introduction of an additional elastic field describing the stress distribution in an equivalent "immobile" packing can explain the observed discrepancy. Work in this direction is currently under progress.

Fig. 1 :

 1 Fig. 1: (Colour on-line) (a)-(c) time and ensemble averaged streamwise velocity vx for a Ω = 6 rpm rotation velocity: (a) map with normalized colormap, (b) with linear and (c) logarithmic z-axis as a function of the depth z at the center of the drum (x = 0). (d) Volume fraction ν as a function of the depth z at the center of the drum (x = 0). Velocities are nondimensionalized by √ gd, where g is the gravity constant and d the mean bead diameter. In (b) and (c), solid lines are the result of a linear and an exponential data fit.

Fig. 2 :

 2 Fig. 2: (Colour on-line) Effective friction coefficient τxz| γd |/ γ d xz P as a function of the inertial number for time and ensemble averaged data, (a) and (b) over the whole 2D drum, (c) and (d) at the center of the drum (x = 0), with linear, (a) and (c), and logarithmic, (b) and (d), x -axis. Green/light gray points correspond to Ω = 12 rpm, red/intermediate gray points to Ω = 6 rpm, and blue/dark gray points to Ω = 2 rpm.

  ) and (b) show the variations of the ratio τ xz | γd |/ γ d xz P as a function of I for time and ensemble averaged data corresponding to the whole 2D drum and for three different rotation velocities. If the rheology were correct, these data would collapse on the single master curve µ(I). This is clearly not the case as we observe an important dispersion of the data -even for the larger values of I. Additionally, as one can see in figs. 2(c) and (d), the collapse predicted by eq. (

  a) |τ | P = µ(I), b) τ ij and γd ij have identical principal directions,

Fig. 3 :

 3 Fig. 3: (Colour on-line) Invariant effective friction coefficient |τ |/P as a function of the inertial number for time and ensemble averaged data, (a) and (b) over the whole 2D drum, (c) and (d) at the center of the drum (x = 0), with linear, (a) and (c), and logarithmic, (b) and (d), x -axis. Same color/gray level code as in fig. 2.

Fig. 4 :

 4 Fig. 4: (Colour on-line) (a) Volume fraction ν and (b) angle difference between principal directions of τij and γd ij as a function of the inertial number for time and ensemble averaged data corresponding to the whole 2D drum. Same color/gray level code as in fig. 2.

Fig. 5 :

 5 Fig. 5: (Colour on-line) (a) Map of the deviation |τ |/P -|τ |/P (x = 0, z) of the effective friction coefficient |τ |/P from the one, µ(I, x = 0), at the center of the drum for Ω = 6 rpm. (b) Angle difference between principal directions of τij and γd ij as a function of the inertial number for Ω = 12 rpm. In (b), data point colors are mapping the x-position in the drum.

Fig. 6 :

 6 Fig. 6: (Colour on-line) (a) Map of the angle difference between principal directions of τij and γd ij for Ω = 12 rpm. (b) Map of the logarithm of the norm of the strain rate tensor times the half drum rotation duration log | γd |/2Ω for Ω = 12 rpm. We indicate with a black line the limit where this deformation δ ≡ | γd |/2Ω = 1.