B D O L L E T

A S C A G L I A R I N I

Two-dimensional plastic flow of foams and emulsions in a channel: experiments and lattice Boltzmann simulations

, have parallels in the behaviour of other soft-glassy materials.

Introduction

Foams and emulsions are dispersions of a fluid phase in a liquid phase, stabilised by surfactants. The dispersed phase is constituted of gas bubbles in foams, and liquid droplets in emulsions. These discrete objects are packed together and jammed, which makes foams and emulsions complex fluids: they exhibit a yield stress σ Y below which they do not flow, but deform elastically. Above yield stress, they flow like rheothinning fluids. Rheometric measurements in a Couette cell or in cone-plate geometry have shown that the shear stress σ and the shear rate γ obey an empirical Herschel-Bulkley law: σ = σ Y + A γn , with A the plastic viscosity and n an exponent generally lower than 1, and often close to 0.5 [START_REF] Princen | Rheology of foams and highly concentrated emulsions. IV. An ex-perimental study of the shear viscosity and yield stress of concentrated emulsions[END_REF][START_REF] Marze | Aqueous foam slip and shear regimes determined by rheometry and multiple light scattering[END_REF][START_REF] Denkov | The role of surfactant type and bubble surface mobility in foam rheology[END_REF], with some dependence on the surfactants used [START_REF] Denkov | The role of surfactant type and bubble surface mobility in foam rheology[END_REF].

The aforementioned measurements did not give access to the micro-structure under flow, and other techniques have been developed to visualise it. In emulsions, confocal microscopy on systems of matched optical index have recently enabled to measure the local structure (Jorjadze, Pontani & Brujić 2013) and the velocity field [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF][START_REF] Mansard | Boundary conditions for soft glassy flows: slippage and surface fluidization[END_REF]. The latter could also be measured using magnetic resonance imaging [START_REF] Ovarlez | Wide-gap Couette flows of dense emulsions: Local concentration measurements, and comparison between macroscopic and local constitutive law measurements through magnetic resonance imaging[END_REF]. In foams index matching is not possible and the route has been to devise bidimensional (2D) experiments on either bubble rafts at the surface of a pool of soap solution, with or without a confining top plate [START_REF] Lauridsen | Velocity profiles in slowly sheared bubble rafts[END_REF][START_REF] Dollet | Two-dimensional flow of foam around an obstacle: Force measurements[END_REF][START_REF] Wang | Impact of boundaries on velocity profiles in bubble rafts[END_REF]Katgert, Möbius & Van Hecke 2008;[START_REF] Katgert | Couette flow of two-dimensional foams[END_REF], or on bubble monolayers confined between two plates in a Hele-Shaw cell (Debrégeas, Tabuteau & di Meglio 2001).

Among many interesting features such as shear banding (see e.g. Schall & Van Hecke (2010) for a review), these studies have called the Herschel-Bulkley law found in rheometry into question. Among possible flow configurations, the Poiseuille flow in a straight channel is particularly interesting, since this geometry enforces a linear variation across the channel of the shear stress, which vanishes at the centre and reaches its maximum at the side walls. Together with an evaluation of the shear rate from the measured velocity profile, it gives access to the relation σ( γ) at the local scale. In particular, [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF] and [START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF] have measured this relation in a series of experiments on emulsions, and they have shown that it did not collapse on a single Herschel-Bulkley law. This deviation from a single flow curve was ascribed to wall effects, more precisely to a nonlocal influence of plastic events happening in the vicinity of the boundaries. The velocity profiles were well fitted by a fluidity model [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF]. This model, based on a kinetic theory approach [START_REF] Bocquet | Kinetic theory of plastic flow in soft glassy materials[END_REF], predicts that the fluidity, defined as f = γ/σ, is proportional to the rate of plastic events and follows a nonlocal diffusion equation when it deviates from its bulk value. The range of influence ξ appearing in this equation, called the spatial cooperativity, was shown to be of the order of a few times (typically, five) the size of the elementary micro-structural constituent (the drop radius in the case of emulsions) [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF][START_REF] Geraud | Confined flows of a polymer microgel[END_REF]. This picture was later applied to other soft materials, such as Carbopol gels [START_REF] Geraud | Confined flows of a polymer microgel[END_REF], granular media [START_REF] Amon | Hot spots in an athermal system[END_REF][START_REF] Kamrin | Nonlocal constitutive relation for steady granular flow[END_REF], and foams in a 2D cylindrical Couette geometry [START_REF] Katgert | Couette flow of two-dimensional foams[END_REF]). The fluidity model agrees with existing experiments, and provides a convenient framework to rationalise the flow of complex fluids. However, at least two points remain unclear and deserve further investigation. The first is the boundary condition at solid walls for fluidity. As a matter of fact, most experimentalists have set it as a free fit parameter, which certainly improves the agreement between the measurements and the predictions from the fluidity model, but does not provide any insight on the role of the walls. Only recently, [START_REF] Mansard | Boundary conditions for soft glassy flows: slippage and surface fluidization[END_REF] explored the role surface boundary conditions for the flow of a dense emulsion. They show that both slippage and wall fluidisation depend non-monotonously on the roughness, a behaviour that has been interpreted with a simple model invoking the building of a stratified layer and the activation of plastic events by the surface roughness. These results are interesting and call for further verification in terms of numerical simulations [START_REF] Benzi R | Rheological properties of soft-glassy flows from hydro-kinetic simulations[END_REF][START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF]) and other complex fluids [START_REF] Katgert | Couette flow of two-dimensional foams[END_REF]. Second, the fluidity parameter f has not been yet convincingly related to an independent measure of the local density of plastic events. In experiments, only indirect indications of such a relation have been proposed, based on the correlations of the fluctuations of the shear rate [START_REF] Jop | Microscale rheology of a soft glassy material close to yielding[END_REF]. Using numerical simulations based on the bubble model [START_REF] Durian | Bubble-scale model of foam mechanics: Melting, nonlinear behavior, and avalanches[END_REF], [START_REF] Mansard | A molecular dynamics study of nonlocal effects in the flow of soft jammed particles[END_REF] were able to measure independently the fluidity and the density of plastic events, but they show that the two quantities are not proportional; more precisely, the rearrangement rate was found to be a sublinear power (with an exponent 0.4) of the fluidity.

Actually, fluidity models offer a potential explanation for the deviation from a unique relation between stress and strain rate, but they are not the only ones. Another approach has been to develop elasto-viscoplastic models (see e.g. [START_REF] Cheddadi | Steady Couette flows of elasto-viscoplastic fluids are non-unique[END_REF] for a review of them) which, in essence, supplement the viscoplastic Herschel-Bulkley rheology by a description of elasticity. These models are local, but since they treat elastic deformation as an independent variable, they also predict deviations from a single Herschel-Bulkley relation. They have been compared with experiments in Couette flows [START_REF] Cheddadi | Steady Couette flows of elasto-viscoplastic fluids are non-unique[END_REF], but not for Poiseuille flows.

All these theoretical approaches rely crucially on the modelling of plastic events, and how they affect the elastic stress and the flow. However, although this connection between elasticity, plasticity and flow has been studied in foam flows in complex geometries [START_REF] Dollet | Two-dimensional flow of foam around a circular obstacle: local measurements of elasticity, plasticity and flow[END_REF][START_REF] Dollet | Local description of the two-dimensional flow of foam through a contraction[END_REF][START_REF] Cheddadi | Understanding and predicting viscous, elastic, plastic flows[END_REF], there is no existing experimental measurement of the rate of plastic events in a Poiseuille flow. 2D foams are particularly well suited for such a study, because elementary plastic events (so-called T1 events) are well characterised by the neighbour swapping of four bubbles (Figs. 2 and3) and are accessible by image analysis, more easily that in other soft-glassy materials.

In this paper, we provide experimental measurements of the rate of plastic events in a Poiseuille flow, on a confined foam in a Hele-Shaw geometry. Such a measurement has already been done in Couette flows (Kabla & Debrégeas 2003;[START_REF] Wang | Impact of boundaries on velocity profiles in bubble rafts[END_REF], but never in a Poiseuille configuration, to our knowledge. We show that it is closely related to the independently measured velocity profiles, and that there is still a non-vanishing plastic activity towards the centre of the channel. The study of the spatial distribution in the number of plastic events and the simultaneous analysis of the velocity profiles allows to bridge between the details of the irreversible plastic rearrangements and the corresponding cooperativity flow behaviour at the macroscopic level [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF][START_REF] Geraud | Confined flows of a polymer microgel[END_REF]. We choose to explore this connection by looking at the relationship between the localisation length of the velocity profiles and the localisation length of the number of plastic events. To complement the cooperativity mechanisms studied in foam with those of "other" soft-glassy systems, we compare the experiments with simulations of emulsion droplets based on the lattice Boltzmann method [START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF]) (LBM). The numerical model falls within the class of "mesoscopic" models, which have been vigorously pursued in the literature to study the behaviour of soft-glassy flows [START_REF] Sollich | Rheology of Soft Glassy Materials[END_REF][START_REF] Sollich | Rheology of Soft Glassy Materials[END_REF]Hébraud & Lequeux 1998;[START_REF] Bocquet | Kinetic theory of plastic flow in soft glassy materials[END_REF][START_REF] Mansard | A molecular dynamics study of nonlocal effects in the flow of soft jammed particles[END_REF][START_REF]Spatial Cooperativity in Microchannel Flows of Soft Jammed Materials: A Mesoscopic Approach[END_REF]. Its elementary, mesoscopic rules, at the level of the lattice units, are designed to reproduce the continuum behaviour of soft-glassy materials, in the very same way that Boltzmann equation in statistical physics, once properly averaged, yields the Navier-Stokes equations. The numerical model we use possesses two advantages that are rarely present together. From one side, it gives a realistic structure of the emulsion droplets, like for example the Surface Evolver method [START_REF] Cox | On the structure of quasi-two-dimensional foams[END_REF][START_REF] Reinelt | Simple shearing flow of dry soap foams with tetrahedrally closepacked structure On the structure of quasi-two-dimensional foams[END_REF][START_REF] Kern | Two-dimensional viscous froth model for foam dynamics[END_REF]; at the same time, due to the built-in properties, the model gives direct access to equilibrium and out-of-equilibrium stresses [START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF], including elastic and the viscous contributions. In contrast to other mesoscopic models, such as Durian's bubble model [START_REF] Durian | Bubble-scale model of foam mechanics: Melting, nonlinear behavior, and avalanches[END_REF], our model naturally incorporates the dissipative mechanisms and the interfacial stresses. Numerical simulations are particularly helpful to go beyond some of the limitations of the experiments: because of wall friction, indeed, there is no simple relation between shear stress and shear rate, whereas the numerical simulations can be performed both with, and without, wall friction. In this paper, wall friction will always refer to the friction of the flowing foam on the confining top and bottom plates (see Fig. 1), and not on the side walls, unless explicitly stated. Numerical simulations also offer the possibility to test the robustness of some of the experimental findings versus a change in the viscous ratio χ between the dispersed phase and the continuous phase, this being set to χ = 1 in all the numerical simulations, whereas χ ≈ 10 -2 in foams; in that sense, the simulations look closer to emulsions. The paper is organised as follows. In Sec. 2, we describe the experimental set-up along with the tools of image analysis for characterisation of the plastic events. In Sec 3, and supplementary material presented in Appendices A and B, we review our computational model based on the LBM. The review of the computational model will be accompanied by further benchmark tests on the capability of the model to include crucial properties as disjoining pressure and wall friction. Results and discussions will be the subject of Sec. 4. The experimentally measured velocity profiles (Sec. 4.1) will be compared with local linear and nonlinear models (Sec. 4.2 and Appendix C). Results of numerical simulations and comparisons with the fluidity model [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Bocquet | Kinetic theory of plastic flow in soft glassy materials[END_REF] will be the subject of Sec. 4.3. In Sec. 4.4, we compare the localisations of the velocity profiles and of the rate of plastic events. In Sec. 4.5, we will finally report details on the orientation of the plastic rearrangements in the flowing material. Conclusions will follow in Sec. 5.

Experimental methods

Setup

We have adapted the setup described in [START_REF] Dollet | Local description of the two-dimensional flow of foam through a contraction[END_REF]. The foam flows in a Hele-Shaw cell, made of two horizontal glass plates of length 170 cm and width 32 cm, separated by a gap h = 2 mm thin enough that the foam is confined as a bubble monolayer (Fig. 1a). Two plastic plates of thickness 2 mm are inserted aside the Hele-Shaw cell, so that the width H of the channel is reduced to 10.66 cm (Fig. 1b). These plates have a negligible roughness compared to the bubble size. The channel is connected upstream to a vertical chamber (Fig. 1a) in which a soap solution is fed at a prescribed flow rate Q l thanks to a syringe pump (PHD2000, Harvard Apparatus). Nitrogen is continuously blown through injectors at the bottom of this chamber, producing rather monodisperse bubbles (Fig. 1b). The flow rate in each injector is independently controlled with an electronic flow-rate controller (Brooks). We identify the liquid fraction ϕ l as the ratio of the liquid flow rate to the total flow rate:

ϕ l = Q l /(Q g +Q l )
, with Q g the gas flow rate. The resulting foam accumulates on top of the chamber, then flows through the channel. The transit time through the whole channel is less than 10 minutes in all experiments; we do not observe significant change of bubble size during this time, hence coarsening is negligible. The soap solution is a mixture of sodium lauryl-dioxyethylene sulfate (SLES), cocoamidopropyl betaine (CAPB) and myristic acid (MAc), following the protocol described in [START_REF] Golemanov | Surfactant mixtures for control of bubble surface mobility in foam studies[END_REF]: we prepare a concentrated solution of 6.6% wt of SLES and 3.4% of CAPB in ultra-pure water, we dissolve 0.4% wt of MAc by continuously stirring and heating at 60 • C for about one hour, and we dilute 20 times in ultra-pure water. The solution has a surface tension Γ = 22.4 mN/m. The contraction region is lit by a circular neon tube, giving an isotropic and nearly homogeneous illumination over a diameter of about 20 cm. Movies of the foam flow are recorded with a CCD camera at a frame rate of 8 frames per second, with an exposure time of 8 ms. The movies are constituted of 1000 images of 1312 × 672 pixels. The pressure drop is measured across the observation zone, by a water-water differential manometer connected to two points of the channel (Fig. 1b) through tubes full of water. We have performed five different experiments, a summary of which parameters is provided in Tab. 1.

Image analysis

To extract the relevant rheological information from the movies, we follow a home-made procedure very similar to that presented in [START_REF] Dollet | Two-dimensional flow of foam around a circular obstacle: local measurements of elasticity, plasticity and flow[END_REF] and [START_REF] Dollet | Local description of the two-dimensional flow of foam through a contraction[END_REF]; we report to these papers for full details. The velocity field is obtained after averaging of all the displacements of all individual bubbles between consecutive frames (about 3 × 10 6 in total). Averaging is performed along 53 lanes aligned with the flow direction. The T1s are tracked as described in [START_REF] Dollet | Two-dimensional flow of foam around a circular obstacle: local measurements of elasticity, plasticity and flow[END_REF]. For the four bubbles concerned by a T1, we denote r d (r a ) the vector linking the centres of the two bubbles that lose (come into) contact, θ ′ d and θ ′ a the angle of these vectors with respect to the flow direction, that we can restrict to the interval [-π/2, π/2] because the orientations of r d and r a are irrelevant (Fig. 2 and Fig. 3), and x d (x a ) the position of the midpoint of the centres of the two bubbles that lose (come into) contact. In our program, the detection of appearing and disappearing contacts is first run independently. As a second step, and3). The quantity dP/dx is the pressure drop along the channel, and dσ/dy is the gradient of shear elastic stress across the channel.

to identify a T1 and minimise artefacts, we decide that a pair of an appearing and a disappearing contact constitute a single T1 if (i) they are on the same image or if the appearing contact happens in the image next to the disappearing one; the latter condition is necessary, because it happens that transient fourfold vertices are erroneously recognised as artificial small bubbles; By following the displacements of the bubbles between subsequent images, we are able to determine the features of a T1 rearrangement. In gray (black) we report the bubble edges just before (after) the T1. With the solid (dashed) line, we report the link between the centres of the two bubbles that lose (come into) contact during the T1. From the analysis of the links, we are able to determine the angles associated with the links that disappear (d) or appear (a) in the T1 rearrangement.

(ii) the positions x d and x a are closer than a critical distance (that we choose to be of the order of the bubble size, to separate from T1s occurring in the neighborhood); (iii) |θ ′ d -θ ′ a | is larger than a critical angle (we choose π/4), this condition being necessary because of the apparition of the aforementioned spurious bubbles. By visual inspection on 30 images, we estimate that this procedure leads to an uncertainty of no more than 5% on the number N T 1 of T1s. We then define the quantity (x d + x a )/2 as the position of a T1, and we ascribe this information to the box where this position belongs. We thus compute the scalar field of the frequency of T1s per unit time and area:

R T 1 = N T 1 2A box t movie ,
where A box is the area of a box and t movie the duration of a movie. Our T1 detection has two major advantages: (i) it is directly based on the topological rearrangements, contrary to indirect characterisations based on velocity correlations; (ii) it yields an unprecedented statistics, up to 2.5 × 10 4 individual T1s, which enables to average over the same lanes as for the velocity and to perform quantitative analysis. We now address the measurement of stress. [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF] has derived the general expression of the stress in a suspension of force-free particles, including the effect of surface tension. This has first been applied to the calculation of the elastic stress in foams and emulsions by [START_REF] Khan | Rheology of foams: I. Theory for dry foams[END_REF], and it yields the same prediction of the elastic shear modulus as [START_REF] Princen | Rheology of foams and highly concentrated emulsions. I. Elastic properties and yield stress of a cylindrical model system[END_REF], which has been experimentally validated (Princen & Kiss 1986) and used ever since. Considering a representative volume element V of foam, the stress writes:

σ = - 1 V ∑ k V k P k 1 -ϕ ℓ p1 + µ V ∫ V ℓ γ dV + Γ V ∫ S (1 -n ⊗ n)dS, (2.1)
where the index k labels the bubbles contained within V , V k and P k are the volume and pressure of bubble k, p is the pressure in the continuous phase of volume V ℓ and viscosity µ, γ is the rateof-strain tensor, S is the surface of the gas/liquid interfaces contained within V and n is the local normal unit vector. The two first terms in (2.1) give the pressure contribution to the stress, and the third term gives the viscous contribution. We cannot measure these contributions in our setup. To systematically analyze plastic events, we perform a Voronoi tessellation from the centres of mass of the droplets. Following the Voronoi tessellation in time, we are able to identify T1 events and associated disappearing (red solid line) and appearing (blue dashed line) links. In gray (black) we indicate the Voronoi cells soon before (after) a T1 event. The numerical results will be compared with the experimental results (see also Fig. 2).

The last term in (2.1) is the elastic stress σ e [START_REF] Cantat | Foams, Structure and Dynamics[END_REF]. A specific advantage of quasi-2D foams, confined in a Hele-Shaw cell so as to form a bubble monolayer, is the possibility to measure elastic stress directly from image analysis. Neglecting the curvature of the films between two bubbles, i.e. assuming that they are flat rectangles of height h and of horizontal side ℓ, which is a reasonable approximation in practice (Fig. 2), the elastic stress equals:

σ e = 2Γhρ ℓ ⟨ ℓ ⊗ ℓ ℓ ⟩ , (2.2)
with ρ ℓ the areal density of bubble edges, and where the average is computed over the same lanes as for the velocity and plastic events. The total number of bubble edges ℓ treated for each experiment is about 10 7 . The "force-free" assumption in the approach of [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF] amounts to neglecting the effect of buoyancy on the bubbles or, more precisely, to neglect the variation of hydrostatic pressure ρgh over the height of the bubbles (with ρ = 10 3 kg/m 3 the density of the soap solution, g = 10 m/s 2 the gravity acceleration, and h = 2 mm the cell gap) compared to the capillary overpressure Γ/R P B , where R P B is the radius of the Plateau borders between three neighbouring bubbles, or two neighbouring bubbles and a plate, the latter ones containing most of the solution for foams confined within a Hele-Shaw cell. The quantity R P B is readily related to the liquid fraction ϕ ℓ : since in order of magnitude, each bubble of volume Ah is surrounded by a total length

≈ √ A of such Plateau borders of cross-section ≈ R 2 P B , one has R P B ≈ h 1/2 A 1/4 ϕ 1/2 ℓ .
Hence, the ratio of the hydrostatic to capillary pressures is

ρghR P B /Γ ≈ ρgh 3/2 A 1/4 ϕ 1/2
ℓ /Γ = 0.01 or 0.02 for our experimental values of the parameters (Tab. 1), much lower than 1. Therefore, it is safe to assume that the bubbles are "force-free" particles.

Numerical method

For the numerical simulations, we adopt a dynamic rheological model based on the lattice Boltzmann method (LBM) [START_REF] Benzi R | The lattice Boltzmann equation: theory and applications[END_REF][START_REF] Doolen | Lattice Boltzmann Method for fluid flows[END_REF][START_REF] Aidun | Lattice Boltzmann Method for complex flows[END_REF]. Historically, the main successful applications of LBM in the context of computational fluid dynamics pertain to the weakly compressible Navier-Stokes equations [START_REF] Benzi R | The lattice Boltzmann equation: theory and applications[END_REF][START_REF] Doolen | Lattice Boltzmann Method for fluid flows[END_REF]) and models associated with more complex flows involving phase transition/separation [START_REF] Chen | Lattice Boltzmann Model for simulating flows with multiple phases and components[END_REF][START_REF] Chen | Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation[END_REF][START_REF] Benzi R | Mesoscopic lattice Boltzmann modeling of soft-glassy systems: Theory and simulations[END_REF]. In particular, we will make use of a computational model for non-ideal binary fluids, which combines a positive surface tension, promoting the formation of diffuse interfaces, with a positive disjoining pressure, inhibiting droplet (or bubble) coalescence. The model has already been validated for a wide spectrum of problems/phenomena [START_REF] Benzi R | Mesoscopic lattice Boltzmann modeling of soft-glassy systems: Theory and simulations[END_REF][START_REF] Benzi R | Herschel-Bulkley rheology from lattice kinetic theory of soft glassy materials[END_REF][START_REF] Benzi R | Rheological properties of soft-glassy flows from hydro-kinetic simulations[END_REF][START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF]). Among the others, these problems/phenomena include the emergence of non-Newtonian Herschel-Bulkley rheology [START_REF] Benzi R | Herschel-Bulkley rheology from lattice kinetic theory of soft glassy materials[END_REF], importance of load-conditions on rheology [START_REF] Benzi R | Rheological properties of soft-glassy flows from hydro-kinetic simulations[END_REF], cooperativity flows [START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF]) and ageing [START_REF] Benzi R | Mesoscopic lattice Boltzmann modeling of soft-glassy systems: Theory and simulations[END_REF]. In this section we just review the method and highlight its essential supramolecular features. The mesoscopic kinetic model considers two fluids A and B, each described by a discrete kinetic distribution function f ζi (r, c i , t), measuring the probability of finding a particle of fluid ζ = A, B at position r and time t, with discrete velocity c i . In other words, the mesoscale particle represents all molecules contained in a unit cell of the lattice. The distribution functions evolve in time under the effect of free-streaming and local two-body collisions, described, for both fluids (ζ = A, B), by a relaxation towards a local equilibrium (f (eq) ζi ) with a characteristic time τ LB :

f ζi (r + c i , c i , t + 1) -f ζi (r, c i , t) = - 1 τ LB ( f ζi -f (eq) ζi ) (r, c i , t) + S (tot) ζi (r, c i , t). (3.1)
Local equilibria are given by a low Mach number expansion of the Maxwellian distribution, namely:

f (eq) ζi = w(|c i | 2 )ρ ζ [ 1 + v • c i c 2 s + vv : (c i c i -c 2 s 1) 2c 4 s ] (3.2)
with w(|c i | 2 ) a set of weights chosen in such a way to maximise the algebraic degree of precision in the computation of the hydrodynamic fields, while c s = 1/ √ 3 is a characteristic velocity (a constant in the model). Our lattice scheme features nine discrete velocities [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF], whose details and associated weights are reported in Tab. 3 in Appendix A. Coarse grained hydrodynamical densities are defined for both species ρ ζ = ∑ i f ζi as well as a global momentum for the whole binary mixture j

= ρv = ∑ ζ,i f ζi c i , with ρ = ∑ ζ ρ ζ .
Non-ideal forces (F ζ ) and a body force term (F b ) are introduced with the source term S (tot) ζi in Eq. (3.1). Such source term produces the desired force once it is projected on the momentum field, i.e.

∑ ζ ∑ i S (tot) ζi c i = ∑ ζ F ζ + F b .
The non-ideal forces include a variety of interparticle forces,

F ζ = F (r) ζ + F (F )
ζ . First, a repulsive (r) force with strength parameter G AB between the two fluids

F (r) ζ (r) = - G AB ρ 2 0 ρ ζ (r) ∑ i=1-8,ζ ′ ̸ =ζ w(|c i | 2 )ρ ζ ′ (r + c i )c i (3.3)
is responsible for phase separation [START_REF] Benzi R | Mesoscopic lattice Boltzmann modeling of soft-glassy systems: Theory and simulations[END_REF]). The parameter ρ 0 is a characteristic normalisation parameter, used as a free parameter in the model. The "short" range interaction in Eq. (3.3) is extended up to energy shells |c i | 2 = 2 (lattice links have been normalised to a characteristic lattice velocity). Phase separating interactions (3.3) are nothing but a lattice transcription of continuum mean-field models for phase segregation [START_REF] Bastea | Hydrodynamics of Binary Fluid Phase Segregation[END_REF]. When the strength parameter G AB /ρ 2 0 in (3.3) is chosen above a critical value, the model achieves phase separation and promotes the emergence of stable diffuse interfaces with a positive surface tension, which allows for the simulation of a collection of droplets. However, the thin films between neighboring droplets are not stable against rupture, as the interactions (3.3) give rise only to negative disjoining pressures [START_REF] Chen | Lattice Boltzmann Model for simulating flows with multiple phases and components[END_REF][START_REF] Chen | Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation[END_REF]. To promote the emergence of a positive disjoining pressure stabilizing the thin films, we introduce a mechanism for frustration (F ) for phase separa-tion with the help of competing interactions [START_REF] Shore | Logaritmic slow domain growth in nonrandomly frustrated systems: Ising Models with competing interactions[END_REF][START_REF] Seul | Domain shapes and patterns: the phenomenology of modulated phases[END_REF]. In particular, we model short range (nearest neighbor, NN) self-attraction, controlled by strength parameters G AA,1 < 0, G BB,1 < 0), and "long-range" (next to nearest neighbor, NNN) self-repulsion, governed by strength parameters G AA,2 > 0, G BB,2 > 0) [START_REF] Chen | Lattice Boltzmann Model for simulating flows with multiple phases and components[END_REF][START_REF] Chen | Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation[END_REF]

F (F ) ζ (r) = -G ζζ,1 ψ ζ (r) ∑ i=1-8 w(|c i | 2 )ψ ζ (r + c i )c i -G ζζ,2 ψ ζ (r) ∑ i=1-24 p(|c i | 2 )ψ ζ (r + c i )c i (3.4) with ψ ζ (r) = ψ ζ [ρ(r)] a suitable pseudo-potential function. The pseudo-potential ψ ζ (ρ ζ ) is taken in the form originally suggested by
ψ ζ [ρ ζ (r)] = ρ 0 [1 -e -ρ ζ (r)/ρ0 ].
(3.5)

The parameter ρ 0 actually marks the density value above which non-ideal effects come into play.

The prefactor ρ 0 in (3.5) is used to ensure that for small densities the pseudopotential is linear in the density ρ ζ . Despite their inherent microscopic simplicity, the above dynamic rules are able to promote a host of non-trivial collective effects [START_REF] Benzi R | Mesoscopic lattice Boltzmann modeling of soft-glassy systems: Theory and simulations[END_REF][START_REF] Benzi R | Herschel-Bulkley rheology from lattice kinetic theory of soft glassy materials[END_REF]. The model gives direct access to the hydrodynamical variables, i.e. density and velocity fields, as well as the local (in time and space) stress tensor in the system, the latter characterised by both the viscous as well as the elastic contributions (see Eq. (A 1) in Appendix A, where both contributions appear). A direct link between the "micro-mechanics" of the model (3.3)-(3.4) and the more familiar "macroscopic" concepts of surface tension and disjoining pressure can actually be established [START_REF] Benzi R | Mesoscopic lattice Boltzmann modeling of soft-glassy systems: Theory and simulations[END_REF][START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF]. To quantify the emergence of the surface tension and the disjoining pressure, one has to consider a 1D problem. For a planar 1D interface, developing along y, the surface tension Γ is a direct consequence of the pressure tensor developing at the non-ideal interface and is computed as the integral of the mismatch between the normal (N) and tangential (T) components of the pressure tensor. Such surface tension scales as as [START_REF] Benzi R | Mesoscopic lattice Boltzmann modeling of soft-glassy systems: Theory and simulations[END_REF])

Γ = ∫ +∞ -∞ [P N -P T (y)] dy ∝ - ∑ ζ=A,B Gζζ ∫ ( dψ ζ dy ) 2 dy - G AB ρ 2 0 ∫ dρ A dy dρ B dy dy. (3.6)
The quantity Gζζ = G ζζ,1 + 12 7 G ζζ,2 comes from a proper combination of the coefficients in the competing interactions. For repulsive interactions, (G AB > 0) the second integral at the right hand side is positive-definite, since (dρ A /dy)(dρ B /dy) < 0. With a proper use of the competing interactions, one can choose Gζζ > 0, and the first term in the right hand side of Eq. (3.6) is negative-definite; consequently, one can decrease the surface tension by simply increasing ρ 0 . The decrease of the surface tension goes together with an increase of the disjoining pressure at the thin film interface. The emergence of a positive disjoining pressure Π d (h) can be controlled in numerical simulations by considering a thin film with two non-ideal flat interfaces, separated by the distance h. Following [START_REF] Bergeron | Forces and structure in thin liquid soap films[END_REF], we write the relation for the corresponding tensions

Γ f (h) = 2Γ + ∫ Π d (h) Π d (h=∞) h dΠ d (3.7)
where Γ f is the overall film tension. Similarly to what we have done for the surface tension Γ, the expression for Γ f is known in terms of the mismatch between the normal and tangential components of the pressure tensor [START_REF] Toshev | Thermodynamic theory of thin liquid films including line tension effects[END_REF][START_REF] Derjaguin | Theory of Stability of Colloids and Thin Films[END_REF],

Γ f = ∫ +∞ -∞ [P N - P T (y)] dy,
where, in our model, P N -P T (y) = p s (y) can be computed analytically [START_REF] Shan | Pressure tensor calculation in a class of nonideal gas lattice Boltzmann models[END_REF][START_REF] Sbragaglia | Interaction pressure tensor for a class of multicomponent lattice Boltzmann models[END_REF]. All the detailed expressions for the interaction stress tensor are reported in Appendix A. From the relation s(h) = Γ f (h) -2Γ it is possible to compute the disjoining pressure: a simple differentiation of s(h) permits to determine the first derivative of the disjoining pressure, ds(h)/dh = h dΠ d /dh. This information, supplemented with the boundary condition Π d (h → ∞) = 0, allows to completely determine the disjoining pressure of the film [START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF]. In Fig. 4 we analyse quantitatively some of these features. In particular we consider the interaction parameters G AB = 0.405, G ζζ,1 = -9.0, G ζζ,2 = 8.1 with ρ 0 chosen in the interval [0.72 : 0.84]. All numbers are reported in lbu (lattice Boltzmann units). As we can see, by increasing the value of ρ 0 , we enhance the energy barrier at the onset of the film rupture.

The body force F b = F P + F D in Eq. (3.1) contains the driving due to the imposed (constant) pressure gradient (F P ) and a drag force (F D ) mimicking the friction between bubbles and confining plastic plates, as in the experimental setup (Fig. 1). Such drag force is taken to be proportional to the velocity vector, as in Janiaud, [START_REF] Janiaud | Two-dimensional foam rheology with viscous drag[END_REF], i.e.

F D = -βv.
(3.8)

Once the droplets are stabilised with a positive disjoining pressure, different packing fractions and polydispersity of the dispersed phase can be achieved. In the numerical simulations presented in the following sections, the fraction of the continuous phase (i.e. the equivalent of the liquid fraction in the foam experiment) is kept approximately equal to ϕ l ≈ 7.5%.

It must be stressed that the micro-mechanics of the model, Eqs. (3.3)-(3.4), is not meant to mimic "specific" physico-chemical details of foams, but rather to model a "generic" soft-glassy model with non-ideal fluid behavior (e.g., nonideal equation of state, phase separation), interfacial phenomena (e.g., surface tension, disjoining pressure) and hydrodynamics. Mesoscale softglassy models are indeed frequently used to characterise the rheology and cooperativity flow of soft-glassy materials [START_REF] Durian | Bubble-scale model of foam mechanics: Melting, nonlinear behavior, and avalanches[END_REF][START_REF] Mansard | A molecular dynamics study of nonlocal effects in the flow of soft jammed particles[END_REF]. As already stressed in the introduction, our numerical model represents a step forward in this direction, in that it provides two basic advantages whose combination is not common. On one hand, it provides a realistic structure of the emulsion droplets, like for instance the Surface Evolver method [START_REF] Cox | On the structure of quasi-two-dimensional foams[END_REF][START_REF] Reinelt | Simple shearing flow of dry soap foams with tetrahedrally closepacked structure On the structure of quasi-two-dimensional foams[END_REF][START_REF] Kern | Two-dimensional viscous froth model for foam dynamics[END_REF]; at the same time, due to its built-in properties, the model gives direct access to dissipative mechanisms in thin films. This latter point will be further discussed and detailed in Appendix B. We also remark that the viscous ratio between the dispersed phase and the continuous phase is kept fixed to χ = 1 (the simulation parameters are summarised in Tab. 2). This choice is dictated by purely numerical reasons, as numerical instabilities emerge when one considers the case of a viscous ratio much smaller or much larger than unity. Nevertheless, we can use this as an advantage in our joint numerical and experimental study, as it offers the possibility to test the robustness of the experimental findings versus a change in the viscous ratio χ between the dispersed phase and the continuous phase. It is also comforting that the latest version of our GPU code [START_REF] Bernaschi | Graphics processing unit implementation of lattice Boltzmann models for flowing soft systems[END_REF]) allows for the simulation of emulsion droplets and their statistics in a reasonable amount of time. The current version runs on multiple-GPU and, by using a combination of CUDA streams and non-blocking MPI primitives, it is able to overlap completely the computation within the bulk of the domain with the exchange of the boundaries. Most simulations have been carried out on Kepler "Titan" GPUs, featuring 14 Streaming Multiprocessors, with a total of 2688 cores running at 0.88 Ghz and a memory bandwidth exceeding 200 GBytes/sec. Each run, spanning multi-million time steps for every single set of parameters, takes less than 12 hours, to be compared with a running time of about 30 hours on previous generation (Fermi) GPU cards. The speedup with respect to a highly tuned (multi-core) CPU version is above one order of magnitude. To develop a systematic analysis of plastic events, we perform a Voronoi tessellation (using the voro++ libraries [START_REF] Rycroft | Analysis of granular flow in a pebble-bed nuclear reactor[END_REF])) constructed from the centres of mass of the droplets, a representation which is particularly well suited to capture and visualise plastic events in the form of droplets rearrangements and topological changes, occurring within the material. Π(h) h F (r) F (r) + F (F) with ρ 0 =0.72 F (r) + F (F) with ρ 0 =0.76 F (r) + F (F) with ρ 0 =0.8 F (r) + F (F) with ρ 0 =0.84

RUN FP |u wall | v0 β β * lbu lbu lbu lbu P1 5 × 10 -7 0.0 2.15 × 10 -2 0 0 P2 5 × 10 -7 0.0 1.05 × 10 -2 10 -5 100 P3 5 × 10 -7 0.0 8.05 × 10 -3 2 × 10 -5 200 P4 4 × 10 -7 0.0 6.10 × 10 -3 2 × 10 -5 200 P5 3 × 10 -7 0.0 4.55 × 10 -3 0 0 P6 3 × 10 -7 0.0 1.65 × 10 -3 2 × 10 -5 200 C1 0.0 2 × 10 -2 0.0 0 0 C2 0.0 2 × 10 -2 0.0 10 -5 100 TABLE 2.

Results and discussions

Experimental velocity profiles

The velocity profiles measured for five experiments are shown in Fig. 5. They are quite flat at the centre (y = 0) of the channel, although not completely flat as would be expected from a Herschel-Bulkley model, and decrease significantly close to the side walls (y = ±H/2). They are well fitted by an exponential profile, which writes either v(y) = v 1 (1 + A cosh y/L v ), or:

v(y) = v 0 cosh(H/2L v ) -α s -(1 -α s ) cosh(y/L v ) cosh(H/2L v ) -1 , (4.1)
with a set of three fitting parameters: v 0 , α s and L v , which have a clear physical meaning: v 0 = v(y = 0) is the centreline velocity, and:

α s = v(y = ±H/2) v 0 (4.2)
is the relative slip, i.e. the ratio of the slip velocity to the centreline velocity. The parameter L v , that we will henceforth call the velocity localisation length, describes the range of influence of the wall friction on the velocity profile. The values of the best fitting parameters are reported in Tab. 1. Among the four experiments run at constant control parameters except the driving flow rate, the relative slip tends to decrease, and the velocity localisation length to increase, at increasing flow rate, except the experiment at flow rate 102.5 mL/min. The fifth experiment is run at larger liquid fraction that the four other: it shows a larger relative slip, and a smaller velocity localisation length, than the experiment with comparable flow rate.

Comparison of experiments with a local model 4.2.1. Linear model

To provide analytical reference equations for the velocity profiles and place our work in the context of the existing literature, we start by comparing our velocity profiles to local models. We start by a comparison to the model of Janiaud, [START_REF] Janiaud | Two-dimensional foam rheology with viscous drag[END_REF], that we adapt to the Poiseuille configuration and to slip boundary conditions. It is appealing, owing to its simplicity, and it has been shown to reproduce well experimental velocity profiles for foam flows in plane Couette geometry (Katgert, Möbius & Van Hecke 2008). The model considers a steady unidimensional flow, where inertia vanishes identically. We also neglect end effects, hence assume that flow is streamwise invariant. Hence, the flow profile writes: v = v(y)e x , with x the streamwise direction and y the spanwise one. The streamwise invariance implies a constant pressure drop: ∇P = e x dP/dx with dP/dx constant. From depth-averaged momentum conservation, 0 = ∇ • σ -∇P + 2f D /h, with f D the foam/wall friction force per unit area. The term 2f D /h is analogous to the drag force F D used in simulations (see Sec. 3). Taking the streamwise component of the equation, we get:

0 = dσ dy - dP dx + 2 h f D , (4.3)
where σ is the xy component of the stress tensor. The model assumes that for the shear stress: σ = σ Y f e (γ/γ Y ) + η γ with γ the shear strain and γ = dv/dy the shear rate, σ Y and γ Y the yield stress and the yield strain, respectively, η a plastic viscosity of the foam, and f e a function quantifying the variation of the elastic stress with the shear strain. Inserting this model in (4.3) yields:

0 = η d 2 v dy 2 + σ Y dγ dy f ′ e ( γ γ Y ) - dP dx -βv,
where f D is assumed to be proportional to the velocity, and β defined as:

f D = - 1 2 hβv. (4.4)
If we neglect the elastic term for simplicity, we get the following ODE for the velocity:

d 2 v dy 2 - v L 2 0 = - v 1 L 2 0 ,
where we have introduced the friction length:

L 0 = √ η β . ( 4.5)
A first boundary condition comes from the fact that x is a symmetry axis, hence v is an even function of y, and we recover the exponential profile (4.1): v(y) = v 1 [1 + A cosh(y/L 0 )], first proposed in Sec. 4.1 as an empirical fit, with the characteristic velocity v 1 proportional to the pressure gradient:

v 1 = - L 2 0 η dP dx . (4.6)
As shown in Sec. 4.1, it turns out that this functional form, with v 1 , A and L 0 as free fitting parameters, reproduces very well the experimental profiles (Fig. 5). However, there is a second boundary condition, coming from a force balance of the foam at the side wall:

σ = ±f D at y = ±H/2, (4.7)
which differs from the original model of Janiaud, [START_REF] Janiaud | Two-dimensional foam rheology with viscous drag[END_REF], who assumed noslip boundary conditions. A macroscopic, visible signature of the balance (4.7) is the angle between the bubble edges and the side walls in Fig. 1b, see also [START_REF] Dollet | Deformation of soap films pushed through tubes at high velocity[END_REF]. Inserting (4.4), (4.5) and (4.6) in (4.3), dσ/dy = β(v -v 1 ), hence:

σ (y = H/2) = βv 1 A ∫ H/2 0 cosh y L 0 dy = βv 1 AL 0 sinh H 2L 0 ,
which we insert in (4.7) to get:

v(y) = v 1 [ 1 - h cosh(y/L 0 ) 2L 0 sinh(H/2L 0 ) + h cosh(H/2L 0 ) ] . (4.8)
This new functional form, with v 0 and L 0 as free fitting parameters, does not fit the experiments. Actually, there is a major discrepancy with the relative slip; setting y = H/2 in (4.8) and using the definition (4.2) of the relative slip, we get:

α s = 1 1 + h 2L0 coth H 2L0 ≃ 1 1 + h/2L 0 , (4.9)
since coth(H/2L 0 ) is very close to 1 for all our experiments. This prediction is much higher than the experimental value, except for the wet foam (Fig. 6).

Nonlinear model

A possible reason for the discrepancy lies in the fact that the wall friction force is nonlinear in velocity, and the bulk viscous stress is nonlinear in shear rate. Following [START_REF] Denkov | The role of surfactant type and bubble surface mobility in foam rheology[END_REF], the bulk viscous stress for a 3D foam equals σ = (τ V F + τ V S )Γ/R (film and interface contribution respectively) with R the bubble radius, and:

τ V F = 1.16Ca 0.47 γ (1 -ϕ l ) 5/6 (0.26 -ϕ l ) 0.1 ϕ 0.5 l , (4.10)
with Ca γ = µ γR/Γ the capillary number (µ = 10 -3 Pa s: bulk viscosity), and τ V S = 9.8πB γ0.18 with B = 2.12 × 10 -3 S.I. an empirical constant for SLES/CAPB/MAc foams. Using as orders of magnitude from the experiments R ≈ √ A/π ≈ 2 mm and γ ≈ v 0 /L v ≈ 1 s -1 , we get τ V S /τ V F ≈ 0.6, hence the film term is dominant, although the surface term is not negligible. Keeping only the film term for simplicity, Eq. (4.10) shows that the bulk viscous stress scales sublinearly with the shear rate:

σ = η ′ γ0.47 , (4.11)
where the prefactor η ′ (primed to distinguish it from the plastic viscosity in the linear law used in Sec. 4.2.1) is:

η ′ = 1.16 µ 0.47 Γ 0.53 R 0.53 (1 -ϕ l ) 5/6 (0.26 -ϕ l ) 0.1 ϕ 0.5 l .
(4.12)

For solutions giving rigid interfaces, like SLES/CAPB/MAc [START_REF] Golemanov | Surfactant mixtures for control of bubble surface mobility in foam studies[END_REF], foam/wall friction is quantified by the force per unit area (or equivalently the wall stress) [START_REF] Denkov | The role of surfactant type and bubble surface mobility in foam rheology[END_REF]:

f D = Γ R [ 1.25C IF √ Ca * √ F 3 1 -F 3 + 2.1C IL (Ca * ) 2/3 ] F 3 , (4.13)
with two empirical constants C IF = 3.7 and C IL = 3.5, Ca * = µv/Γ another capillary number, and:

F 3 = √ 1 -3.2 ( 1 -ϕ l ϕ l + 7.7
) -1/2 .

For v ≈ 1 cm/s, the ratio of the second term to the first term in (4.13) is 5, hence we neglect the second term. Eq. (4.13) then shows that the wall friction force scales sublinearly with the velocity:

f D = - 1 2 hβ ′ √ v, (4.14)
with the following value of the wall friction constant (primed to distinguish it from its counterpart in the linear law used in Sec. 4.2.1):

β ′ = 2.5C IF hR √ Γµ √ F 3 3 1 -F 3 . (4.15)
Like in Sec. 4.2.1, see Eq. ( 4.5), we can construct a characteristic length from η ′ and β ′ . To do so, it is convenient to replace the exponent 0.47 by 1/2 in (4.11), recasting the factor γ0.03 in the definition (4.12) of η ′ ; this factor is almost constant, and equal to 1, for all our experiments. The characteristic length is then

L ′ 0 = (η ′ /β ′ ) 2/3
. It is the extension the friction length L 0 defined in (4.5), to the case of the nonlinear wall friction (4.14) and the nonlinear bulk viscous stress (4.11). We compute with all the experimental values of the parameters appearing in (4.12) and (4.15): L ′ 0 = 1.9 mm for ϕ l = 4.8%, and 2.7 mm for ϕ l = 16.9%. These orders of magnitude are compatible with the experimental values of the localisation lengths (Tab. 1).

The effect of nonlinear wall friction and bulk viscous stress on the velocity profile has been theoretically considered for Couette flows [START_REF] Weaire | Velocity dependence of shear localisation in a 2D foam Phil[END_REF][START_REF] Weaire | A simple analytical theory of localisation in 2D foam rheology Phil[END_REF], but not for Poiseuille flows. Therefore, in Appendix C, we compute analytically the velocity profile using the nonlinear laws (4.11) and (4.14), and we show that the role of these nonlinearities on the velocity is negligible.

Hence, the present model is too simple to capture wall slip in our experiments, which suggests that the role of elastic stresses is crucial. This is qualitatively supported by the fact that the only experiment for which the local model is quite accurate in predicting the amount of slip is for a wet foam, which stores less elastic energy [START_REF] Cantat | Foams, Structure and Dynamics[END_REF]. To further support this idea, we plot the shear component of the elastic stress and the normal elastic stress difference in Fig. 7; see the end of Sec. 2.2 for the measurement of the elastic stress. The shear elastic stress is indeed about four times weaker for the wet foam than for the four other experiments. For these experiments at given bubble area and liquid fraction, its variation across the channel is as follows: towards the centre of the channel, although with a significant asymmetry for some experiments, there is a zone of quasilinear increase around σ e xy = 0. The width of this region decreases slightly at increasing flow rate. Outside this region, the elastic shear stress plateaus to a value which does not depend much on the flow rate. Interestingly, there is still some velocity variation, and a significant plastic activity, outside those regions where the shear elastic stress plateaus. Except the experiment for the wet foam, the normal elastic stress difference σ e xx -σ e yy is always positive, i.e. the bubbles are elongated streamwise, an effect which is clearly visible on Fig. 1b. It tends to increase towards the wall.

Eq. (4.3) expresses the balance between the driving pressure gradient, the foam/wall friction, and the gradient of elastic and bulk viscous stresses. Close to the middle of the channel, the velocity gradient is very weak, hence bulk viscous stress is negligible, and the gradient of the shear elastic stress is roughly constant (Fig. 7 Tab. 1; the pressure gradient is always larger than the gradient of shear elastic stress, the missing part being wall friction. This is a major difference with Poiseuille experiments in 3D channels [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF][START_REF] Geraud | Confined flows of a polymer microgel[END_REF], where wall friction is absent. This prevents us from measuring directly the spanwise stress from the pressure gradient, contrary to the aforementioned studies.

Numerical simulations and comparison with the fluidity model

We have shown the inaccuracies of a local model without elasticity to capture our experimental data thanks to the inspection of the boundary condition at the wall. To test the effects of elasticity, local visco-elastoplastic models could be used [START_REF] Cheddadi | Steady Couette flows of elasto-viscoplastic fluids are non-unique[END_REF]), but it is not straightforward to deduce from them testable predictions. Moreover, elasticity is not the only aspect of the physics of foams and dense emulsions missed by a model such as that discussed in the previous sections. It has been recently shown in experiments of flowing emulsions in microchannels [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF]) that in order to capture the velocity profiles, the intrinsic non-local rheology of such soft-glassy materials must be considered. Non-locality is due to the long-ranged relaxation of stress released after a plastic rearrangement occurs. Building up on previous models for soft-glassy rheology [START_REF] Sollich | Rheology of Soft Glassy Materials[END_REF][START_REF] Sollich | Rheology of Soft Glassy Materials[END_REF]Hébraud & Lequeux 1998), [START_REF] Bocquet | Kinetic theory of plastic flow in soft glassy materials[END_REF] proposed a kinetic elasto-plastic model (KEP henceforth) describing explicitely spatial interactions among plastic events. The continuum limit of KEP suggests that the local relation between stress σ and rate of strain γ is dictated by a kind of inverse effective viscosity, the fluidity

f = γ σ , (4.16)
which obeys a non-local diffusion-relaxation equation of the form

ξ 2 ∆f (x) + f b (σ(x)) -f (x) = 0, (4.17)
where f b is the bulk fluidity, i.e. the value of the fluidity in the absence of spatial cooperativity (ξ = 0) and ξ is the cooperativity length scale within the material due to spatial heterogeneities and represents, basically, a correlation length for the fluidity itself, as it can be easily derived from Eq. (4.17). The other important result of KEP is that such fluidity must be expected to be proportional to the rate of occurrence of plastic events R T 1 in the system, i.e. f ∝ R T 1 . Furthermore, KEP encompasses the effect of elasticity through the nonlocal relaxation of elastic stress induced by plastic events. There is, however, a difficulty in testing this nonlocal model against our Tab. 2); the experimental velocity profile with flow rate 152.5 ml/min (see Tab. 1) is also reported to show that we are able to tune the wall friction parameters in the numerics to achieve the same localisation observed in the experiments. The equivalent β * in the experiments is β * ≈ 250 (see text for details). On the abscissae, the y-location across the channel has been normalised by the total channel height.

experiments. The role of wall friction is crucial in experiments, whereas the nonlocal model has been set up and tested in its absence, although recent studies have considered the coupled role of wall friction and nonlocality [START_REF] Barry | Nonlocal effects in the continuum theory of shear localisation in 2d foams Phil[END_REF][START_REF] Scagliarini | Non-locality and viscous drag effects on the shear localisation in soft-glassy materials[END_REF]. Indeed, wall friction complicates the stress profile across the channel, as discussed in Sec. 4.2, and it is thus not straightforward to extract relevant flow curves σ( γ) from our experiments. Hence, it is interesting to run numerical simulations, where the wall friction can be set off and tuned at will. Various sets of numerical simulations have been performed in the (F P , β * ) parameter space (see Tab. 2 for the numerical values used), where β * is meant to be the value of friction coefficient (Eq. 3.8) β made dimensionless with the channel width H and viscosity η, i.e.

β * = βH 2 η = H 2 L 2 0 , (4.18)
where the last equality is based on the definition of the friction length L 0 given in (4.5). A flat velocity profile in the bulk is shown by all curves (Fig. 8), including the case with β * = 0, witnessing the presence of a non trivial bulk rheology [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF][START_REF] Geraud | Confined flows of a polymer microgel[END_REF]. We also report the experimental velocity profile with flow-rate 152.5 ml/min (see Tab. 1), just to show that we are able to tune the wall friction parameters in the numerics to achieve the same localisation observed in the experiments for which an equivalent wall friction parameter β * ≈ 250 can be estimated based on the wall friction constant (4.15) and the plastic viscosity (4.12). While a direct comparison of experiments and numerics in terms of velocity profiles is complicated by boundary conditions (as observed in the experimental data of Fig. 5, slippage is found to occur at the surfaces of the experiments, while the numerical simulations are performed by imposing no-slip at the walls), an important insight will be provided by simulations in validating the picture of the plastic flow (Sec. 4.4) at changing the wall friction constant which is a freely tunable parameter in the numerical model, unlike in experiments.

The top panel of Fig. 9 indeed provides some indications that wall friction does not seem to dramatically affect the distribution of plastic events. There we plot the rate of plastic rearrangements, normalised by the total number of events, from experiments and numerics (for the three β * 's). Data show a moderately good collapse onto each other. At a given driving pressure drop, increasing wall friction results in a decrease of the total number of plastic events N T1 . We could estimate the number N T1 in the numerical simulations and it is reported in the bottom panel of Fig. 9: for the same simulation time (see caption of Tab. 2) plastic events diminish from N T1 ∼ 6 × 10 3 to N T1 ∼ 2 × 10 3 for increasing β * from 0 to 400. A similar trend is observed for the centreline velocity which is reported in the inset of the bottom panel of Fig. 9. To make this statement more quantitative, we notice that the overall decrease in the number of plastic events can be well captured by the function

G(β * ) = 2b 2 N (0) T 1 β * [ 1 - 1 cosh( √ β * /b) ] , (4.19)
a scaling behaviour that can be obtained from the expression of the centreline velocity, v 0 in (4.1), with α s = 0 (no-slip boundary condition for the numerics). The parameter

N (0)
T 1 in equation (4.19) sets the number of plastic events in the limit β * → 0 whereas the argument √ β * /b of the hyperbolic function in Eq. (4.19) is inversely proportional to the velocity localisation length. The choice of equation ( 4.19) as a fitting function is suggested by the consideration that the total number of plastic events is dominated by events occurring in boundary regions where the shear stress is approximately constant and, hence, the fluidity f is basically proportional to the shear rate γ. Consequently, being the number of events, by definition, equal to the integral of the corresponding rate R T 1 and since

R T 1 ∝ f ≈ | γ|, we can assume that G(β * ) ∝ ∫ H/2 0
| γ|dy, which equals the centreline velocity. Interestingly, the estimate of b that we get from a best fit procedure (b ≈ 6.0) is greater than the estimate of b based on the friction length L 0 in (4.5), which would yield b = 2. This is an indication that the velocity localisation length in the numerical simulations is larger than the localisation length induced by the wall friction F D = -βv. This is not a surprise, because our numerical simulations had already confirmed the presence of a cooperativity length scale [START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF], without wall friction, in a Couette flow setup. This supports the idea that an effective localisation length results from the sum of the friction length plus the cooperativity length [START_REF] Barry | Nonlocal effects in the continuum theory of shear localisation in 2d foams Phil[END_REF], an issue that we will further explore in Sec. 4.4.

The study of the spatial distribution in the number of plastic events and the simultaneous analysis of the localisation in the velocity profiles, allows to bridge between the "microscopic" details of local irreversible plastic rearrangements and the macroscopic flow. A connection between the rate of T1 events and the fluidity field is indeed visible in the top panel of Fig. 9. The dashed line indicates sinh(y/L v )/y, which is the "synthetic" fluidity profile (up to an unessential numerical scaling factor) based on the hyperbolic cosine fit of the velocity profile (see Sec. 4.1) and a linear variation of the shear stress across the channel. Interestingly, a significant plastic activity remains towards the centre of the channel, and it is well correlated to the fluidity field, which remains finite in such regions, whereas the strain-rate goes to zero. Moreover, a closer inspection reveals that the decrease in the number of plastic events is affected by the wall friction constant β * , with a steeper decrease associated with the larger β * . Fig. 9 calls therefore for a deeper understanding with regard to the link between the rate of plastic events and the local flow properties.

Localisation lengths: comparison of plasticity and shear rate

To go further, we choose to explore the connection between the rate of plastic events and the local flow properties, by looking at the relationship between the localisation length of the velocity profiles, L v , and the localisation length of the number of plastic events, L p (henceforth called plastic localisation length). This connection enables to compare experiments and simulations, despite their different boundary conditions. The velocity localisation length L v is estimated by a hyperbolic cosine function cosh(y/L v ), from which the decay length L v is extracted (see Sec. 4.1). Simultaneously, the plastic localisation length L p is computed out of an exponential fit of the symmetrised rate of plastic events close to the wall (Top Panel of Fig. 10). Since our numerical simulations have already confirmed the presence of a cooperativity length scale [START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF] without wall friction, they are good candidates to complement the experimental findings, showing how the spatial distribution of plastic events is affected by a change in the wall friction parameter β. Hence, in the Bottom Panel of Fig. 10 we also look at the localisation in the numerics, by fixing the pressure gradient and changing β * , something that cannot be easily done in experiments with the data at hand. At fixed pressure gradient, we show the log-lin plot of the rate of plastic events from simulations with different β * . The extracted L p is found to be a decreasing function of β * . 1. These data have been symmetrised with respect to the centreline; the error bars denote the standard deviation resulting from this averaging. Data below 10 -4 mm -2 •s -1 are not shown because they are statistically irrelevant (less than ten T1s counted per bin per experiment). Dashed lines represent the linear fits of the data. Bottom Panel: data from numerical simulations complement the experimental results reported in the Top Panel. In particular, to appreciate the effect of wall friction at fixed pressure gradient, we show the log-lin plot of the rate of plastic events from simulations P1-3 (fixed pressure drop and different β * ) close to the bottom wall (the dashed lines represent best linear fits of the data). Inset: Plastic localisation length as function of the wall friction parameter β * .

In Fig. 11 we report a scatter plot of the velocity localisation length L v versus the plastic localisation length L p for three sets of data: experiments (symbols as in Tab. 1), simulations with fixed pressure drop and various β * 's (filled squares) and simulations with fixed β * = 200 and various pressure drops (filled circles). Fig. 11 shows that the two localisation lengths are indeed close to each other. The fact that the values L p and L v agree, confirms the picture of the "plastic flow"; it is also compatible with the fact that the rate of plastic events and the fluidity seem to be proportional (Sec. 4.3). [START_REF] Barry | Nonlocal effects in the continuum theory of shear localisation in 2d foams Phil[END_REF] have combined the local model presented in Sec. 4.2 with the nonlocal equation for the fluidity field (4.17), in the case of a Couette flow with linear laws for the bulk viscous stress and wall friction (as in Sec. 4.2.1). They predicted that the velocity localisation length is an increasing function of both the cooperativity length ξ, and of the friction length L 0 defined by (4.5). Here, this theoretical prediction can be tested for the first time versus our experiments and simulations. Some care is required in doing so because of the nonlinear nature of wall friction and bulk viscous stress in experiments and of the difference between Couette and Poiseuille flows. However, we have shown in Appendix C that the effect of nonlinear wall friction and bulk viscous stress is very weak. Since the velocity localisation length L v is much smaller than the channel width and the stress does not vary much across the localisation zone, the comparison with the Couette predictions is relevant. For this reason, we also repeated some numerical simulations in a Couette flow geometry (see Tab. 2). The associated data nicely collapse on the same master curve, stressing even more the robustness of our findings at changing the load conditions.

In experiments, at given liquid fraction and bubble area, the localisation length L v increases with increasing flow rate. Moreover, Tab. 1 shows that at given flow rate (up to 5%) and bubble area, the localisation length is lower for a wet foam than for a dry one, in qualitative agreement with the experiments of [START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF] on emulsions. We tried to get further insight on these findings, by investigating the local effect of single plastic events, namely by measuring the way they affect displacement and elastic stress fields on their surroundings, in the spirit of [START_REF] Picard | Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids[END_REF] in theory and [START_REF] Desmond | Topological rearrangements and stress fluctuations in quasi-two-dimensional hopper flow of emulsions[END_REF] in experiments. Our hope was to directly evidence the cooperativity length as the range of influence on displacement and stress fields of single plastic events, and whether it would depend on flow rate and liquid fraction. However, the data turned out to be too noisy to address this specific question, in particular because of the difficulty to properly substract the effect of the mean flow. Hence, we cannot directly measure the cooperativity length. In some sense, the work presented here bypasses the problem of an accurate measurement of the cooperativity length, but directly explores the link between localisation phenomena in the velocity profiles and the rate of plastic events.

Orientation of the plastic events

The importance of plastic rearrangements has been stressed in that the occurrence of these events induces long range correlations within the soft-glassy material. It is also acknowledged FIGURE 12. Normalised distributions of the orientations of plastic events. Distributions of appearing (left Panel) and disappearing (right Panel) centre-to-centre links of bubbles involved in rearrangements are shown, from experiment (thin blus line) and simulations (thick red line). θ denotes the angle formed by the link and the direction of the flow, i.e. the positive x-axis. Following [START_REF] Princen | Rheology of foams and highly concentrated emulsions. I. Elastic properties and yield stress of a cylindrical model system[END_REF], the extreme values are found for a dry foam at liquid fraction ϕ l = 0 and are indicated with a dashed line. [START_REF] Picard | Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids[END_REF][START_REF] Schall | Structural rearrangements that govern flow in colloidal glasses[END_REF]) that T1s possess a non trivial angular structure with a quadrupolar topology. It seems, then, reasonable to argue that for a full understanding of the way they determine nonlocal effects inside the system, not only the distribution of their locations in space, but also their orientational properties need to be addressed. Therefore, we go further with the description of plastic events, and study their angular statistics from experiments and simulations. More precisely, focusing on the four bubbles involved in a T1, we define as a disappearing link the segment connecting the centres of the two bubbles which were in contact before the event (and which are then far apart), and as an appearing link the connector between the other two bubble centres (see also Sec. 2 and Figs. 2 and3); we then measure for each event the angle between such links and the flow direction. We have observed that the angles are reversed between both sides of the channel, consistently with the fact that y = 0 is an axis of symmetry. Therefore, we choose to analyse the statistics of the quantities θ d = θ ′ d sign (y) and θ a = θ ′ a sign (y) (see Fig. 1 for the sign convention of y). We did not observe a significant variation of the distribution of these angles across the channel, hence we analyse the distributions of these angles for all T1s, whatever their location across the channel. Fig. 12 shows the histogram of θ d and θ a for one experiment and one simulation, while the average and standard deviation of these quantities are summarised for all experiments in Tab. 1. This analysis shows that T1s have preferential orientations: θ d is peaked around 0.5 rad, with a small dispersion, and θ a around -0.7 rad, with a larger dispersion. The average values do not depend significantly on the flow rate. For the wet foam, θ d is larger, and θ a slightly smaller.

We now derive some reference values for these angles from a micro-structural analysis. Since our foams are rather monodisperse, it is interesting to use the simple geometrical model of a sheared 2D hexagonal foam [START_REF] Princen | Rheology of foams and highly concentrated emulsions. I. Elastic properties and yield stress of a cylindrical model system[END_REF]) (see also [START_REF] Khan | Rheology of foams: I. Theory for dry foams[END_REF]). In this model, the unit cell drawn in dashed lines in Fig. 13 is sheared, and the location of the vertices is computed to comply with the equilibrium rule that the three edges meet at equal angles. To account for the finite liquid fraction, the vertices are decorated with Plateau borders which radius R P is an increasing function of the liquid fraction (Fig. 13, left): ϕ l = (2 √ 3 -π)R 2 P /A h , with A h the area of one hexagon. This structure can be sheared up to the point where two neighbouring Plateau borders meet, which defines the onset of the T1 (Fig. 13,right).

The two angles θ ′ a and θ ′ d can be computed from simple geometry, when the two Plateau borders come into contact (Fig. 13, right). The length of the edge c between the two bubbles about to detach is equal to R P . Now at a given strain γ, this length equals [START_REF] Khan | Rheology of foams: I. Theory for dry foams[END_REF]:

c = a(1 -γ √ 3/2)/ √ 4 + γ 2
, where a is the side length of the undeformed hexagon. Setting c = R P in the latter equation yields the strain γ c at which the T1 occurs; γ c is a decreasing function of R P . Moreover, γ = 1/ √ 3-cot α [START_REF] Princen | Rheology of foams and highly concentrated emulsions. I. Elastic properties and yield stress of a cylindrical model system[END_REF]), and we compute from geometrical considerations in the right panel of Fig. 13:

cot θ ′ d = 2/ √ 3 -cot α = 1/ √ 3 + γ c , and cot θ ′ a = -2/ √ 3 -cot α = - √ 3 + γ c .
Qualitatively, these two expressions show that both θ ′ d and θ ′ a are decreasing functions of γ c , hence increasing functions of R P , hence of the liquid fraction. The extreme values are found for a dry foam at ϕ l = 0, for which γ c = 2/ √ 3 [START_REF] Princen | Rheology of foams and highly concentrated emulsions. I. Elastic properties and yield stress of a cylindrical model system[END_REF], and for the jamming transition for which γ c = 0: θ ′ d varies between π/6 ≃ 0.52 rad (dry foam) and π/3 ≃ 1.05 rad (jamming transition), and θ ′ a between -π/3 and -π/6. Our measured values are indeed in these ranges. The values of the disappearing angles for the four experiments with ϕ l = 4.8% are compatible (within experimental dispersion) with the dry foam prediction, the latter indicated with a vertical dashed line in Fig. 12. The predicted increase of the angles with liquid fraction is compatible with the experiments for θ d , but not for θ a .

Although the model by [START_REF] Princen | Rheology of foams and highly concentrated emulsions. I. Elastic properties and yield stress of a cylindrical model system[END_REF] gives useful reference values, it is difficult to make a more quantitative comparison based on liquid fraction, because the distribution of liquid is specific to each system. In simulations, the films between droplets are thick, and contain a significant proportion of the liquid. In experiments, the distribution of water is complex because of the 3D structure of the bubbles; there is relatively more water close to the confining plates than in the midplane in between [START_REF] Cox | On the structure of quasi-two-dimensional foams[END_REF]. The hexagonal foam model of [START_REF] Princen | Rheology of foams and highly concentrated emulsions. I. Elastic properties and yield stress of a cylindrical model system[END_REF] is a good approximation of the structure of our experimental foams across the midplane between the top and bottom confining plates, but the liquid fraction across this plane, relevant in the hexagonal model, is significantly lower than the experimental liquid fraction. Moreover, the measurement of the appearing angle is less precise than that of the disappearing one, because the relaxation of the four bubbles after a T1 is fast; hence, the measurements made on the image after the topological rearrangement may not be representative of the configuration at the instant of a T1. This also explains why the dispersion is larger for θ a than for θ d .

Conclusions

We have reported on the first experimental study measuring the rate of plastic events in Poiseuille flows of foams. Experiments have been supplemented by numerical simulations, capable to capture the realistic foam structure and to incorporate naturally the expected mesoscopic dynamics. We have addressed the relation between T1 distribution and macroscopic rheology and revealed a link between the localisation lengthscale of the velocity profiles and that of plastic events across the channel, confirming the relevance of cooperativity for foams [START_REF] Katgert | Couette flow of two-dimensional foams[END_REF]). The use of numerical simulations allowed to study in a controlled way (something not easily feasible in the experiments) the effect of wall friction, helping to confirm its role in the emergence of an extra localisation for the velocity profiles, as predicted theoretically [START_REF] Barry | Nonlocal effects in the continuum theory of shear localisation in 2d foams Phil[END_REF]. Our study highlighted that the elasticity gives rise to a complex near-to-wall dynam-ics which calls for focused studies both experimentally (in the spirit of the recent work by [START_REF] Mansard | Boundary conditions for soft glassy flows: slippage and surface fluidization[END_REF]) and numerically, and for a more refined theoretical modelling of the boundary conditions. Finally, unprecedented results on the distribution of the orientation of plastic events show -with good agreement between experiments and numericsthat there is a non-trivial correlation with the underlying local shear strain; this suggests that more complex forms for the propagators invoked in theoretical models of soft-glassy materials [START_REF] Bocquet | Kinetic theory of plastic flow in soft glassy materials[END_REF]) may be needed, with an explicit angular structure, especially in situation of non-homogeneous stress (as it is for Poiseuille flows).

The authors kindly acknowledge funding from the European Research Council under the EU Seventh Framework Programme (FP7/2007-2013) / ERC Grant Agreement no [279004]. We acknowledge computational support from CINECA (IT). MS and AS gratefully acknowledge M. Bernaschi for computational support and F. Bonaccorso for helpful visualisations of the flowing emulsions from the numerical simulations.

Appendix A. Pressure tensor in LBM simulations

In this Appendix we provide the technical details for the lattice Boltzmann pressure tensor used in equations (3.6) and (3.7) to compute both the surface tension and the disjoining pressure at the non-ideal interface. Given the mechanical model for the lattice interactions described in (3.3)-(3.5), an exact lattice theory is available [START_REF] Shan | Pressure tensor calculation in a class of nonideal gas lattice Boltzmann models[END_REF][START_REF] Sbragaglia | Interaction pressure tensor for a class of multicomponent lattice Boltzmann models[END_REF] which allows to connect the interaction forces to the lattice Pressure Tensor. The exact pressure tensor is given by

P αβ = ∑ ζ,i f ζi c α i c β i + ∑ ζ P (int) ζ,αβ . ( A 1)
The term ∑ ζ,i f ζi c α i c β i represents an internal contribution to the pressure tensor and its outof-equilibrium contribution gives the dissipative stress in the system [START_REF] Gross | Shear stress in nonideal fluid lattice Boltzmann simulations[END_REF], while P (A 2)

The contribution coming from the phase separating interactions P

(r) ζ,αβ is (Sbragaglia & Belardinelli 2013) P (r) ζ,αβ = G AB 2 ρ ζ (r) ∑ i=1-8 w(|c i | 2 )ρ ζ ′ (r + c i )c α i c β i ζ ′ ̸ = ζ (A 3)
while the contributions coming from the frustrating interactions are given by various terms,

P (F,1) ζ,αβ , P (F,2) ζ,αβ , P (F,4) ζ,αβ , P (F,5) ζ,αβ , P (F,8)
ζ,αβ , labeled with the number of the "energy shell" (see Tab. 3)

P (F,1) ζ,αβ = G ζζ,1 2 ψ ζ (r) ∑ i=1-4 w(1)ψ ζ (r+c i )c α i c β i + G ζζ,2 2 ψ ζ (r) ∑ i=1-4 p(1)ψ ζ (r+c i )c α i c β i (A 4) P (F,2) ζ,αβ = G ζζ,1 2 ψ ζ (r) ∑ i=5-8 w(2)ψ ζ (r+c i )c α i c β i + G ζζ,2 2 ψ ζ (r) ∑ i=5-8 p(2)ψ ζ (r+c i )c α i c β i (A 5) P (F,4) ζ,αβ = G ζζ,2 4 ∑ i=9-12 w(4)ψ ζ (r)ψ ζ (r+c i )c α i c β i + G ζζ,2 4 ∑ i=9-12 w(4)ψ ζ ( r - c i 2 ) ψ ζ ( r + c i 2 ) c α i c β i (A 6)
) ψ ζ ( r + c i 2 ) c α i c β i .
2009). We start by presenting benchmark computations for the motion of droplets in confined channels at changing the capillary number. The drag force on a single bubble that slides past a solid wall was first investigated by [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] and has recently received renewed attention [START_REF] Denkov | Foam-wall friction: Effect of air volume fraction for tangentially immobile bubble surface Colloids Surf[END_REF][START_REF] Denkov | Viscous Friction in Foams and Concentrated Emulsions under Steady Shear[END_REF][START_REF] Katgert | Flow in linearly sheared twodimensional foams: From bubble to bulk scale[END_REF]. For a single bubble sliding past a solid wall, [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] showed that the drag force scales non-linearly with the capillary number, Ca, defined in terms of the dynamic viscosity of the carrier liquid and the relative velocity between the bubble and the wall. The lubrication approximation yields the velocity U of the bubble immersed in the Poiseuille flow to be approximated at leading order in U/U av = 1 + 1.29(3Ca) 2/3 , (B 1)

where U av represents the mean flow velocity. This theory can be readily modified for the 2D bubble placed in the Poiseuille flow between two parallel plates (Afkhami, Leshansky & Renardy 2011) U/U av = 1 + 0.643(3Ca) 2/3 . (B 2)

An extension to droplets with an arbitrary viscosity has been considered in various papers [START_REF] Schwartz | On the motion of bubbles in capillary tubes[END_REF][START_REF] Hodges | The motion of a viscous drop through a cylindrical tube[END_REF]. For a very viscous drop, results analogous to the previous equations can readily be found, as the coefficients 1.29 and 0.643 in equations (B 1) and (B 2), respectively, are reduced by a factor of 2 -1/3 ≈ 0.794, yielding U/U av = 1 + 1.023(3Ca) 2/3 (B 3) for a bubble in a cylindrical capillary and U/U av = 1 + 0.511(3Ca) 2/3 (B 4) for a 2D bubble in a channel. In Fig. 14 we present our benchmark tests for equation (B 4). As we can see, the predicted scaling for the velocity of the bubble well agrees with the theoretical prediction, confirming a scaling exponent in the capillary number close to 2/3 and a numerical coefficient in between the case of a very viscous droplet and the case of a bubble We now continue by presenting benchmark tests for the drag force between two bubbles sliding past each other, F bb . Some recent works [START_REF] Denkov | Foam-wall friction: Effect of air volume fraction for tangentially immobile bubble surface Colloids Surf[END_REF][START_REF] Denkov | Viscous Friction in Foams and Concentrated Emulsions under Steady Shear[END_REF][START_REF] Katgert | Flow in linearly sheared twodimensional foams: From bubble to bulk scale[END_REF]) have provided evidences that the viscous drag force scales like F bb ∝ Ca ξ , with a scaling exponent ξ between 1/2 and 2/3. This is an important test for our numerical simulations: LBM modelling of two phase flows is intrinsically a diffuse interface method and involves a finite thickness of the interface between the two liquids and related model parameters. The values of the interface thickness and capillary number need to be larger than the one suggested by physical considerations in order to make the simulations affordable [START_REF] Komrakova | Lattice Boltzmann simulations of drop deformation and breakup[END_REF][START_REF] Magaletti | The sharp-interface limit of the Cahn-Hilliard/Navier-Stokes model for binary fluids[END_REF]. Nevertheless, the structure and the dynamical properties of the emulsion droplets that we reproduce in the numerical simulations share nontrivial features with the experiments [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Goyon | How does a soft glassy material flow: finite size effects, nonlocal rheology, and flow cooperativity Soft Matter[END_REF][START_REF] Mansard | Boundary conditions for soft glassy flows: slippage and surface fluidization[END_REF]. It is therefore of great importance to investigate the scaling laws associated with friction properties, to show that they are realistic and in line with those measured in experiments. In particular, we measure the viscous drag forces between bubbles directly by rheological experiments where two rows of ordered bubbles are sheared past each other. Results are reported in Fig. 15. A scaling law in the velocity difference U between the two rows of droplets is confirmed, F bb ∼ Ca ξ , with a scaling exponent between 1/2 and 2/3.

Appendix C. Velocity profile in the nonlinear local model

In this Appendix, we provide a solution of the velocity profile, obeying momentum conservation (4.3) with the nonlinear law (4.14) for wall friction, and (4.11) for the bulk viscous stress FIGURE 14. Velocity of a 2D droplet in a confined channel. The viscous ratio between the dispersed phase and the continuous phase is set to χ = 1 in all the numerical simulations. In the top panel we report two snapshots associated with two different capillary numbers. Blue/white (dark/light) colors indicate regions with majority of the dispersed/continuous phase. The droplet is driven by a constant pressure gradient. The average velocity of the droplet is normalised with respect to the mean flow velocity (Uav) in the inlet of the channel. The scaling laws for both very viscous (χ ≫ 1) droplet [START_REF] Schwartz | On the motion of bubbles in capillary tubes[END_REF][START_REF] Hodges | The motion of a viscous drop through a cylindrical tube[END_REF]) and a bubble in a 2D channel are reported [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF]Afkhami, Leshansky & Renardy 2011). The velocity scaling well agrees with the theoretical prediction, confirming a scaling exponent in the capillary number close to 2/3 and a numerical coefficient in between the two extreme cases (χ ≫ 1 and χ ≪ 1).

Φ l =6% Φ l =10% Φ l =14% Φ l =18%
Φ l =10%, t=15 000 (lbu) Φ l =10%, t=30 000 (lbu) Φ l =10%, t=45 000 (lbu) Φ l =18%, t=15 000 (lbu) Φ l =18%, t=30 000 (lbu) Φ l =18%, t=45 000 (lbu) but with a velocity exponent 1/2 instead of 0.47, an approximation which enables to provide an analytical solution. Hence we solve for -H/2 ≤ y ≤ 0 (so that dv/dy ≥ 0):

0 = η ′ d dy [ ( dv dy ) 1/2 ] - dP dx -β ′ v 1/2 ,
with boundary conditions: dv/dy = 0 at y = 0, and σ = -f D at y = -H/2. We make space and velocity dimensionless: ȳ = y/L ′ 0 and v = v/V , by rescaling them by

L ′ 0 = (η ′ /β ′ ) 2/3 V = [(-dP/dy)/β] 2 . (C 1)

FIGURE 1 .

 1 FIGURE 1. (a) Sketch of the side view of the setup. (b) Snapshot of an experiment. The distance between the two side walls is H = 106.6 mm. The average bubble size is 12.3 mm 2 and the liquid fraction is ϕ l = 4.8%. The two spots at the left and right of the image are the points between which the pressure drop is measured.

FIGURE 2 .

 2 FIGURE 2. Left Panel: A snapshot of the bubbles in the foam flowing from left to right. Right Panel: Sketch of a plastic event.By following the displacements of the bubbles between subsequent images, we are able to determine the features of a T1 rearrangement. In gray (black) we report the bubble edges just before (after) the T1. With the solid (dashed) line, we report the link between the centres of the two bubbles that lose (come into) contact during the T1. From the analysis of the links, we are able to determine the angles associated with the links that disappear (d) or appear (a) in the T1 rearrangement.

FIGURE 3 .

 3 FIGURE 3. Left Panel: A snapshot of the droplets (identified by their corresponding Voronoi cells) in a concentrated emulsion, flowing from left to right, obtained in numerical simulations based on the lattice Boltzmann models. Right Panel: sketch of a T1 plastic event from the simulations. To systematically analyze plastic events, we perform a Voronoi tessellation from the centres of mass of the droplets. Following the Voronoi tessellation in time, we are able to identify T1 events and associated disappearing (red solid line) and appearing (blue dashed line) links. In gray (black) we indicate the Voronoi cells soon before (after) a T1 event. The numerical results will be compared with the experimental results (see also Fig.2).

  Summary of the simulations parameters. The first six rows refer to runs in the Poiseuille (P#) flow setup, while the last two are relative to the Couette (C#) flow numerical simulations. Other relevant parameters (kept fixed among the various runs) are the fraction of the continuous phase ϕ l ≈ 7.5% and the viscous ratio between the dispersed and continuous phase χ = 1. The interaction parameters for the phase separating interactions (see Eq. (3.3)) and competing interactions (see Eqs. (3.4)) are given in the text and the pseudo-potential reference density is ρ0 = 0.83. The disjoining pressure for these interaction parameters is characterised in Fig.4. The total integration time is Ttot = 2 × 10 6 lbu (lattice Boltzmann units) of which Tss = 1.25 × 10 6 steps in the steady state.

FIGURE 4 .

 4 FIGURE 4. This figure shows the emergence of the disjoining pressure in the lattice Boltzmann model (see Eq. (3.1) with phase separating interactions obtained with a repulsive (r) force (see Eq. (3.3)) supplemented with competing interactions (see Eq. (3.4)) whose role is to provide a mechanism for frustration (F ). The use of phase separating interaction is associated with a negative disjoining pressure. Competing interactions stabilise thin films with the emergence of a positive disjoining pressure, the latter tunable with the parameter ρ0 in Eqs. (3.3) and (3.4). Further details can be found in[START_REF] Sbragaglia | The emergence of supramolecular forces from lattice kinetic models of non-ideal fluids: applications to the rheology of soft glassy materials[END_REF].

FIGURE 5 .

 5 FIGURE 5. Velocity profiles for the five experiments described in Tab. 1. These data have been symmetrised with respect to the centreline; the error bars denote the standard deviation resulting from this averaging. Each dashed line is a fit of the data by the law (4.1); see Tab. 1 for the values of the best fitting parameters. The dotted line is a fit of the data series △ with the law (C 4) accounting for nonlinear wall friction and bulk viscous stress; see Sec. 4.2.2 and Appendix C for details. It is barely distinguishable from the dashed line.

FIGURE 6 .

 6 FIGURE 6. Relative slip measured in the experiments described in Tab. 1, as a function of the relative slip (4.9) predicted by the local model.

FIGURE 7 .

 7 FIGURE 7. Left panel: shear component of the elastic stress, and right panel: normal difference of the elastic stress, for the five experiments described in Tab. 1.

FIGURE 8 .

 8 FIGURE 8. Numerical velocity profiles, normalised by the centreline velocity v0, at changing the wall friction constant. Data from three sets of simulations are shown with dimensionless wall friction parameter (see equation 4.18) β * = 0, 100, 200 (runs P1-3 in table Tab. 2); the experimental velocity profile with flow rate 152.5 ml/min (see Tab. 1) is also reported to show that we are able to tune the wall friction parameters in the numerics to achieve the same localisation observed in the experiments. The equivalent β * in the experiments is β * ≈ 250 (see text for details). On the abscissae, the y-location across the channel has been normalised by the total channel height.

FIGURE 9 .

 9 FIGURE 9. Top Panel: plot of the rate of plastic rearrangements as a function of y: experiments (△ in Tab. 1) are compared with numerical data from runs P1-3. The dashed line indicates the function sinh(y/Lv)/y, representing the fluidity profile based on the hyperbolic cosine fit of the velocity profile (see Sec. 4.1 for details). Numerical data have been symmetrised. Bottom Panel: Total number of T1s as function of β * from the simulations; the dashed line is a fit with the functional form given in Eq. (4.19). Inset: the centreline velocity vs β * is reported.

FIGURE 10 .

 10 FIGURE 10. Top Panel: Decimal logarithm of the density of T1s per unit time and area as a function of y/H, for the five experiments described in Tab. 1. These data have been symmetrised with respect to the centreline; the error bars denote the standard deviation resulting from this averaging. Data below 10 -4 mm -2 •s -1 are not shown because they are statistically irrelevant (less than ten T1s counted per bin per experiment). Dashed lines represent the linear fits of the data. Bottom Panel: data from numerical simulations complement the experimental results reported in the Top Panel. In particular, to appreciate the effect of wall friction at fixed pressure gradient, we show the log-lin plot of the rate of plastic events from simulations P1-3 (fixed pressure drop and different β * ) close to the bottom wall (the dashed lines represent best linear fits of the data). Inset: Plastic localisation length as function of the wall friction parameter β * .

FIGURE 11 .

 11 FIGURE 11. Scatter plot of the velocity localisation length Lv (computed from a hyperbolic cosine fit of the velocity profiles) vs the plastic localisation length Lp (computed out of an exponential fit of the symmetrised rate of plastic events across the channel) for three sets of data: experiments (symbols as in Tab. 1), simulations of Poiseuille flow with fixed pressure drop and various normalized friction coefficients β * 's (filled squares) and with fixed β * = 200 and various pressure drops (filled circles) and simulations of Couette flow at two β * 's (filled triangles); both lengths are normalised by the mean bubble diameter. The dashed line is the Lv = Lp curve.

FIGURE 13 .

 13 FIGURE 13. Portion of an unsheared (left) and sheared (right) hexagonal foam. The strain is defined as γ = 4∆x/3a.

  is a contribution coming from the interactions and embeds the elastic stresses of the system. As for P (int)ζ,αβ , we can separately write the contributions coming from the repulsive (r) phase separating interactions (see equation (3.3)), and those coming from competing interactions providing a mechanism of frustration (F) (see equation (3

  4). The first belt lattice velocities are indicated with i = 1...8 while the second belt ones with i = 9...24. p(|ci| 2 ) or w(|ci| 2 ) indicate the weight associated with the i-th link in the various interactions. The weights associated to the velocity at rest, w(0) and p(0), are chosen to enforce a unitary normalisation,[ψ ζ (r + c 5 )ψ ζ (r + c 3 ) + ψ ζ (r + c 1 )ψ ζ (r + c 7 )] c α 13 c β 13 + + G ζζ,2 4 p(5) [ψ ζ (r + c 5 )ψ ζ (r + c 4 ) + ψ ζ (r + c 2 )ψ ζ (r + c 7 )] c α 14 c β 14 + + G ζζ,2 4 p(5) [ψ ζ (r + c 2 )ψ ζ (r + c 8 ) + ψ ζ (r + c 6 )ψ ζ (r + c 4 )] c α 15 c [ψ ζ (r + c 6 )ψ ζ (r + c 1 ) + ψ ζ (r + c 3 )ψ ζ (r + c 8 )]

FIGURE 15 .

 15 FIGURE 15. Viscous drag force F bb between droplets measured directly in rheological experiments where two rows of ordered bubbles are sheared past each other. The packing fraction of the continuous phase into the dispersed phase is changed in the interval ϕ l = [0.06 : 0.18]. The left panel reports three snapshots of the simulations for two different packing fractions. Blue/white (dark/light) colors indicate regions with majority of the dispersed/continuous phase. A scaling law in the Capillary number is found, F bb ∼ Ca ξ , with a scaling exponent between 1/2 and 2/3.

TABLE 1 .

 1 Summary of the main characteristics of the experiments presented in this paper, with the respective symbols used in the figures. The parameters v0, αs and Lv come from the fit of the velocity profiles by the formula (4.1). The angles θ d and θa are the orientations of plastic events(Figs. 2

	flow rate ϕ l	bubble	v0	Lv	αs	θ d	θa	dP/dx dσxy/dy symbol
	(ml/min) (%) area (mm 2 ) (mm/s) (mm)		(rad)	(rad)	(kPa/m) (kPa/m)
	27.5	4.8 13.0 ± 2.4	2.3	3.3 0.62 0.54 ± 0.14 -0.79 ± 0.24 1.08	0.50	•
	52.5	4.8 12.6 ± 2.3	4.3	3.9 0.43 0.46 ± 0.16 -0.67 ± 0.33 1.08	0.68
	102.5	4.8 12.1 ± 2.1	8.5	8.7 0.34 0.46 ± 0.19 -0.75 ± 0.46 2.23	0.51	♢
	152.5	4.8 12.3 ± 1.7 12.4	6.8 0.31 0.46 ± 0.17 -0.73 ± 0.41 2.27	0.91	△
	160.2 16.9 15.1 ± 1.9 12.9	4.6 0.90 0.73 ± 0.22 -0.86 ± 0.21 0.91	0.087	▽

TABLE 3

 3

	Phase Separating Interactions Eq. (3.3) Shell	lattice Links
	w(|ci| 2 )		|ci| 2	ci
	4/9		0	(0, 0)
	1/9		1 (±1, 0); (0, ±1)
	1/36		2	(±1, ±1)
	1st term, rhs Eq. (3.4) 2nd term, rhs Eq. (3.4) Shell	lattice links
	w(|ci| 2 )	p(|ci| 2 )	|ci| 2	ci
	4/9	247/420	0	(0, 0)
	1/9	4/63	1	(±1, 0); (0, ±1)
	1/36	4/135	2	(±1, ±1)
	0	1/180	4	(±2, 0); (0, ±2)
	0	2/945	5 (±2, ±1); (±1, ±2)
	0	1/15120	8	(±2, ±2)

. Links and weights of the two belts, 25-speeds lattice

[START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF][START_REF] Benzi R | Mesoscopic lattice Boltzmann modeling of soft-glassy systems: Theory and simulations[END_REF]

) for all interactions given in equations (3.3) and (3.

Appendix B. Friction forces in LBM simulations

In this Appendix, we propose benchmark tests for the LBM introduced in section 3 with regard to the viscous drag forces acting on individual bubbles. Friction properties in the thin films between neighbouring droplets/bubbles or between droplets/bubbles and the walls are important for both foams and concentrated emulsions [START_REF] Denkov | Foam-wall friction: Effect of air volume fraction for tangentially immobile bubble surface Colloids Surf[END_REF][START_REF] Denkov | Viscous Friction in Foams and Concentrated Emulsions under Steady Shear[END_REF][START_REF] Denkov | The role of surfactant type and bubble surface mobility in foam rheology[END_REF]Katgert et al. Then the problem becomes:

where the subscript designs derivation, with boundary conditions: vȳ = 0 at ȳ = 0, and v1/2

The velocity thus obeys an autonomous equation of the form:

√ y. This kind of equation can be recast as a firstorder ODE [START_REF] Polyanin | Handbook of Exact Solutions for Ordinary Differential Equations[END_REF] by setting w = vȳ : it then becomes wv

The latter is an ODE with separable variables, which is thus simply integrated to yield: 2 w3/2 /3 = -2v + 4v 3/2 /3 + const. The boundary condition at ȳ = 0 imposes that w = 0 for the unknown centreline velocity v0 , hence 2 w3/2 /3 = -2(v -v0 )+4(v 3/2 -v3/2 0 )/3. In the limit 2L ′ 0 /H ≪ 1, which is a good approximation in our experiments, the dimensionless bulk viscous stress term (v

3) The right-hand side is a decreasing function over [0, 1], equal to 1 for v = 0 and to 0 for v = 1. Therefore, for a given value of the parameter h/2L ′ 0 , the boundary condition at the wall: vȳ = h 2 v/4L ′2 0 , admits a single solution for v and vȳ at ȳ = -H/2L ′ 0 . The velocity field obeying (C 2) and the boundary condition then obeys the implicit equation:

)] ,

where 2 F 1 designs the hypergeometric function. With this complicated expression, fitting the experimental data is not easy. A simpler alternative consists in developing (C 3) for 1 -v ≪ 1. This gives: vȳ = (3/4) 2/3 (1 -v) 4/3 , with general solution: 1 -v = 48/(ȳ + const.) 3 . This gives an alternative fitting formula for the velocity profile:

Taking V , L ′ 0 and y 0 yields a fit which is close to the exponential fit (4.1), the difference between both fits being within the dispersion of the experimental data. This suggests that the effect of the nonlinearities of wall friction and of the bulk viscous stress on the flow profile is weak. Finally, we have checked that the qualitative conclusions brought by Fig. 6, namely that the relative slip is overestimated by the local model, remains valid with nonlinear laws. Hence, the role of nonlinearities is secondary in this study, and can be neglected.