H C H R A I B I

R W U N E N B U R G E R

J P E T I T 1 A N D J-P. D E L V I L L E

Eddies and interface deformations induced by optical streaming

Keywords: interfacial flows, electrohydrodynamics effects, boundary integral method

We study flows and interface deformations produced by the scattering of a laser beam propagating through non absorbing turbid fluids. Light scattering produces a force density resulting from the transfer of linear momentum from the laser to the scatterers. The flow induced in the direction of the beam propagation, called "optical streaming", is also able to deform the interface separating the two liquid phases and to produce wide humps. The viscous flow taking place in these two liquid layers is solved analytically, in one of the two liquid layers with a stream function formulation, as well as numerically in both fluids using a Boundary Integral Element Method. Quantitative comparisons are shown between the numerical and analytical flow patterns. Moreover, we present predictive simulations dedicated to the effects of the geometry, of the scattering strength and of the viscosities, on both the flow pattern and the deformation of the interface. Theoretical arguments are finally put forth to explain the robustness of the emergence of secondary flows in a two-layer fluid system.

Introduction

Streaming flow induced by light scattering has been emphasized by [START_REF] Savchenko | Transfer of momentum and torque from a light beam to a liquid[END_REF] in various liquids such as liquid crystals and colloids. It was showed that a light beam having an inhomogeneous transverse intensity profile may enforce convection in liquids as a consequence of the transfer of linear and angular momentum from light to liquids. In a recent investigation on the effect of a continuous laser wave on soft liquidliquid interfaces [START_REF] Schroll | Liquid transport due to light scattering[END_REF], we proposed an original experiment to optically induce bulk flows in turbid liquids without any heating. We showed that light scattering by refractive index fluctuations can produce a steady flow due to transfer of momentum from light to the liquid. In addition, the viscous stress exerted by this flow can deform the liquid interface. We successfully compared Figure 1. Eddies produced by the scattering of a continuous laser beam propagating in a non absorbing turbid liquid (micellar phase of microemulsion) layer bounded by a top glass wall and a bottom transparent brine solution; a thin layer of metastable foam is present at the brine-microemulsion interface. A beam centered toroidal eddy, evidenced by the advection of foam in this picture, induces a centered deformation of the interface by viscous withdrawal in the direction of beam propagation. The interface deformation is represented by the black curve. The experiment is performed in a so-called Winsor II equilibrium [START_REF] Kellay | Properties of surfactant monolayers in relation to microemulsion phase behaviour[END_REF]) obtained here by mixing equal volumes of n-Heptane and brine with a small amount of AOT surfactant (B. P. [START_REF] Binks | Stability of oil-in-water emulsions in a low interfacial tension system[END_REF]. Sodium chloride is used to screen electrostatic repulsion between surfactant heads and thus to significantly reduce the interfacial tension. For the chosen concentrations, [N aCl] = 0.05M and [AOT ] = 40mM , the interfacial tension is σ = 1 -5 10 -6 N/m. The beam power is P = 1.32W and the beam waist ω0 = 2.99µm. Eddies are magnified by overlaying 100 frames from a video sequence.

experimentally observed deformations of the two-fluid interface to their analytical prediction. We call this flow hereafter "optical streaming" by analogy with the acoustic streaming phenomenon [START_REF] Nyborg | Acoustic streaming near a boundary[END_REF]) that occurs when an acoustic wave propagates through a sound absorbing liquid. Indeed, while acoustic streaming is dissipative in nature due to wave absorption, "optical streaming" results from the scattering of the incident beam on density inhomogeneities. However, the net result in terms of induced flow is very similar. More recently, we demonstrated that optical streaming is responsible for the droplet emission at the tip of liquid jets triggered by radiation pressure [START_REF] Wunenburger | Fluid flows driven by light scattering[END_REF]. Optical streaming can be observed experimentally as illustrated in Figure 1. In the present investigation, we consider theoretically the optical streaming taking place in a set of two horizontal liquid layers separated by an interface. They are submitted to a laser beam propagating upward along the z-axis (see figure 2). Assuming that at least one of the layers is a turbid liquid, i.e. contains refractive index inhomogeneities smaller that the optical wavelength, a bulk force is induced in this layer by light scattering leading to the formation of steady toroidal eddies. If the interfacial tension is sufficiently weak, these eddies are able to deform the liquid-liquid interface by the induced viscous stress. The goal of the present investigation is thus to study optical streaming in a more general frame, by analyzing the influence of the physical properties of the fluids and the geometry of the container. Beyond understanding the nature of these flows, an interesting prospect is to use them in future microfluidic experiments [START_REF] Delville | Laser microfluidics: fluid actuation by light[END_REF]).

Our article is structured as follows: the physical modeling of the problem is presented in section 2. In section 3, we analytically model the optical streaming in a single liquid layer, using a stream function formulation. The numerical simulation in a two layer system is presented in section 4, using a Boundary Element Method. In section 5, comparisons between numerical and analytical predictions are discussed, addressing the influence of the aspect ratio and sizes of the cell containing the fluids, of the scattering intensity and of the viscosity ratio on the flow patterns. Finally, good quantitative agreements between numerical simulations and experiments are reported in the appendix confirming the validity of the physical description and the numerical solution proposed in this work.

Formulation of the problem

Geometry and liquid properties

We consider two liquid layers labeled by i, i=1,2, separated by an interface, illuminated by a continuous laser beam, of Gaussian intensity profile along the radial coordinate and of cylindrical symmetry along the vertical propagation axis z. The two liquids are enclosed in a cylindrical container of solid boundaries S W 1 and S W 2 , of radius R ≃ 10 -3 m >> ω 0 ≃ 10 -5 m and thicknesses L 2 ≃ L 1 ≃ 10 -4 m >> ω 0 , where ω 0 is the characteristic laser beam radius, also called beam waist (see figure 2). Considering the axisymmetry of the system along the z axis, we choose cylindrical coordinates (e r , e Θ , e z ) with origin O located at the intersection of the beam axis with the initially flat interface located therefore at z = 0. We refer to a point x by the space coordinates (r, Θ, z). Physical properties of liquids i = 1, 2 are their viscosities η i , densities ρ i , refractive indices n i and momentum attenuation coefficients ̟ i which represent the relative loss of light linear momentum per unit propagation length. The interfacial tension is denoted by σ. The characteristic values of the liquids are

η i ≃ 10 -3 Pa s, ρ i ≃10 3 kg m -3 , n i ≃ 1.33, ̟ i ≃ 70m -1 and σ ≃ 10 -6 N/m.

Bulk force associated to light scattering

The momentum flux associated to light of celerity c in vacuum, propagating in a fluid of refractive index n i is n i c I. The beam intensity profile is

I(r, z) ≈ I(r) = 2P πω 2 0 e -2(r/ω 0 ) 2 (2.1)
P is the beam power.

When light is elastically scattered by the refractive index inhomogeneities of the medium, the forward component of its momentum density q iz decreases along the direction of propagation. To quantify this decrease, we introduce the forward momentum attenuation coefficient in liquid i defined by

̟ i = - 1 q iz ∂q iz ∂z , i = 1, 2 (2.2)
We assumed statistical isotropy of the distribution of the refractive index inhomogeneities within the medium. Therefore, the total rate at which the momentum density decreases due to scattering is ̟ i n i c Ie z , i = 1, 2. As a consequence of momentum conservation, the liquids are subjected to a body force density F i such as

F i = ̟ i n i c Ie z = F i e -2(r/ω 0 ) 2 e z , i = 1, 2 (2.3) 2.3. Flow equations
Since the characteristic order of magnitude of optically driven flows as measured in the experimental investigation reported in [START_REF] Wunenburger | Fluid flows driven by light scattering[END_REF] is u i ≃10 -4 m s -1 , we investigate a viscous flow at low Reynolds number

(Re i = u i L i ρ i /η i ≃ 10 -2 ).
Both liquids satisfy the mass conservation and Stokes equations,

∇ • u i = 0 ; i = 1, 2 (2.4) 0 = -∇p i + η i ∆u i + F i ; i = 1, 2
(2.5) The gravitational force has been included in the corrected pressure term p i that is defined by p i = p ′ i + ρ i gz, p ′ i being the fluid pressure. Assuming no slip at the interface and immiscible fluids, the hydrodynamic stress balance on the interface (described by its height h(r)) involving interfacial tension and gravity effects is written as

B.C.1 : T 1 • n -T 2 • n = (σκ -(ρ 1 -ρ 2 )gh)n on S I (2.6)
Here, T i = -p i I + η i (∇u i + ∇u i t ) is the corrected hydrodynamic stress tensor, n is the unit vector normal to the interface directed from fluid 1 to fluid 2 and

κ(r) = 1 r d dr r dh dr 1 + dh dr 2
is the double mean curvature of the axisymmetric interface in cylindrical coordinates. The continuity of velocity at the interface yields

B.C.2 : u 1 = u 2 on S I (2.7)
At the container walls, the classical no-slip condition for viscous flow is assumed, i.e.

B.C.3

:

u 1 = 0 on S W 1 ; B.C.4 : u 2 = 0 on S W 2
(2.8) Finally, the motion of the interface is described using a Lagrangian approach dx dt = u(x) on S I (2.9)

Analytical solution of a one-fluid model

When the top fluid is much less viscous and much less dense than the bottom one, a simplified analytical solution to the above problem can be derived. Indeed, in this situation of an inviscid top fluid remaining at a constant uniform pressure, the stress exerted on the bottom layer can be neglected leading to a one-fluid model. To develop the solution, the flow taking place in the bottom fluid layer can be solved analytically using a stream function formulation. We assume that the width of the laser beam is much smaller than the width of the resulting eddy (in figure 1 for instance, the eddies presented are ten times larger than the beam waist). Therefore, the scattering force can be represented as a line force f 1 e z = F 1 πω 2 0 e z along the centre line of the cylindrical container. The amplitude of the deformation is supposed to be very small compared to its characteristic width and a flat interface is assumed.

The flow within the layer is modeled as a steady viscous flow contained in a cylindrical layer of height L 1 and radius R centered on the beam axis. Due to the symmetry of the optical excitation, the flow is axisymmetric, and can be described by a stream function ψ(r, z) such that

u = ∇ × ψ(r, z) r e Θ (3.1)
This choice for the stream function ensures incompressibility (Equation (2.4)), and the Stokes equation becomes

0 = -∇p + η 1 ∇ 2 ∇ × ψ(r, z) r e Θ + F 1 (3.2)
Taking the curl of Equation (3.2) removes the pressure gradient, and yields

η 1 D 4 ψ (r, z) =r (∇×F 1 ) • e Θ (3.3)
where

D 4 • ≡ D 2 (D 2 •) and D 2 ≡ ∂ 2 ∂r 2 -1 r ∂ ∂r + ∂ ∂z 2 .
Assuming a line force, the right-hand term of Equation (3.3) is a delta-function at r = 0. Thus, except at r = 0

D 4 ψ (r, z) = 0 (3.4)
Assuming separation of variables, the general solution is

ψ (r, z) = G(r)sin kz (3.5)
where k = π/L 1 . Note that this solution always satisfies ∂ur ∂z (r, z = 0) = 0 and ∂ur ∂z (r, z = -L 1 ) = 0. Even though not strictly correct, this last condition ensures a straightforward analytical solution. The solution to the problem is finally achieved using four boundary conditions. We require by symmetry u r (0, z) = 0, no flux through the far wall of the container r = R, u r (R, z) = 0 and no slip along it u z (R, z) = 0. Replacing Equation (3.5) into Equation (3.4) we find

G (r) = - U f 4πk r K 1 (kr) kr kR I 1 (x) K 1 (x) xdx -I 1 (kr) kr kR K 1 (x) 2 xdx - U f 4πk rA(k) I 1 (kr) kr kR I 1 (x) K 1 (x) xdx -K 1 (kr) kr kR I 1 (x) 2 xdx (3.6)
where

A(k) = kR 0 I 1 (x)K 1 (x)xdx kR 0 I 1 (x) 2 xdx
, I 1 and K 1 being the first order modified Bessel functions of first and second kind respectively. The last free parameter U f can now be estimated by matching the flow amplitude and the forcing. This means that the integrated delta-function forcing must be related to the scattering induced flow. More specifically, we take the cross-sectional integral of Equation (3.3) at the mid-plane of the liquid layer, i.e. at z = -L 1 /2. On the left-hand side, we get

η 1 D 4 ψdS = η 1 U f (3.7)
On the right-hand side, combining Equations (2.1) and (2.3), we get

r (∇×F 1 ) • e Θ dS = 2π ∞ 0 r 2 - ∂F 1 (r) ∂r dr = 2P ̟ 1 n 1 c = f 1 (3.8)
Thus, the characteristic velocity of the light-induced bulk flow U f is

U f = f 1 η 1 = F 1 πω 2 0 η 1 (3.9)
The velocity components are deduced from Equations (3.1) and (3.5) It yields

u r (r, z) = - G(r) r kcos kz (3.10) u z (r, z) = - U f 4π -K 0 (kr) kr kR I 1 (x) K 1 (x) xdx -I 0 (kr) kr kR K 1 (x) 2 xdx sin kz - U f 4π A(k) I 0 (kr) kr kR I 1 (x) K 1 (x) xdx + K 0 (kr)
p(r, z) = η 1 U f 4π k [K 0 (kr) + A(k)I 0 (kr)] cos kz (3.12)
where I 0 and K 0 are the 0 th order modified Bessel functions of first and second kind respectively. The solution of this simplified model allows to express the corrected normal stress at the top of the fluid layer which is given by the T zz component of the corrected stress tensor

T zz (r, z) = -p (r, z) + 2η 1 ∂u z ∂z (3.13)
Using Equations (3.11) and (3.12), we can write

T zz (r, z) = -η 1 U f π kS(r)cos kz (3.14) with S (r) = K 0 (kr) 4 1 -2 kr kR I 1 (x) K 1 (x) xdx - I 0 (kr) 2 kr kR K 1 (x) 2 dx + A(k) 4 I 0 (kr) 1 + 2 kr kR I 1 (x) K 1 (x) xdx + 2K 0 (kr) kr kR I 1 (x) 2 xdx (3.15)
To predict the deformation, the corrected normal stress exerted by the viscous flow in the lower layer is balanced with buoyancy and capillarity. Since the interface deformation is considered to remain small (i.e. h ′ = dh/dr << 1), the corrected hydrodynamic stress experienced by the interface at radius r is well approximated by -T zz (r, 0)e z . The linearized equilibrium equation of the interface is therefore written in a dimensionless form using L 1 as a characteristic length and σ/L 1 as a characteristic pressure

h′′ (r) + h′ (r) r -Bo h(r) + f 1 σ S(r) = 0 (3.16)
where Bo = ρ 1 gL 2 1 /σ is the gravitational Bond number and the symbol . refers to a dimensionless quantity. This equation is solved using a differential algebraic equation solver (DASSL) based on a backward differentiation formula in the interval [r = ǫ, r = R]. Note that ǫ = 0 since T zz (0, 0) diverges due to the line force distribution exerted at r = 0. Since Equation (3.16) is of second order, we require two boundary conditions for h(r). We therefore assume h(r = ǫ) = h ǫ and need the condition dh/dr(r = ǫ) = h ′ ǫ , h ǫ being our shooting parameter. The convergence condition is chosen to be dh/dr(r = R) = 0 with a tolerance of 0.1%. Analytical profiles predictions are compared with experimental and numerical results in the appendix.

To check the predictibility of the present simplified one-layer analytical model, we also investigate numerically the two-fluid configuration. This is presented in the following section.

Numerical resolution

The Boundary Element Method reveals to be very accurate for solving interfacial flow problems with high resolution as reported in the analysis of flow involving electric and magnetic fields [START_REF] Sherwood | Breakup of fluid droplets in electric and magnetic fields[END_REF], optical radiation pressure (Chraibi et al. 2008b[START_REF] Chraibi | Optohydrodynamics of soft fluid interfaces: Optical and viscous nonlinear effects[END_REF] or buoyancy [START_REF] Manga | Low reynolds number motion of bubbles, drops and rigid spheres through fluid-fluid interfaces[END_REF][START_REF] Koch | Numerical and theoretical solutions for a drop spreading below a free fluid surface[END_REF]. When used to solve problems involving surface forces like capillarity, radiation pressure, etc.. and volume forces, which can be written under a conservative form, this method allows to reduce the formulation of the problem to boundary integrals only. In the present case, the scattering force obviously does not belong to this class, and consequently the integral formulation of the problem involves a volume integral term corresponding to the contribution of this force [START_REF] Occhialini | Boundary intergral/spectral element approaches to the navier-stokes equations[END_REF][START_REF] Issenmann | Unsteady deformations of a free liquid surface caused by radiation pressure[END_REF]. The integral formulation is achieved by first making use of the Green's function for the Stokes operator which corresponds to the solution of the fundamental Stokes problem [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF]. The integral equations are presented below under a dimensionless form by choosing U σ = σ/η 1 as a reference velocity, L 1 as a reference length and p i0 = η i σ/(L 1 η 1 ) i = 1, 2 as reference pressures in both fluids. Taking into account the boundary conditions on the interface S I and on the walls S W 1 , S W 2 (see figure 2), and denoting by n e the outward unit normal vector to these boundaries, the integral form providing the velocity at any point x in the interior of the volumes V 1 and V 2 , is given by ūi

(x) = S I +S W i U • T i • n e dS y - S I +S W i n e • K • ūi dS y + V i U • F i dV y i = 1, 2 (4.1)
The solution of this equation requires a prior determination of stress on S I , S W 1 and S W 2 as well as of the velocity on S I . This is achieved by solving an integral equation which, once B.C.1 (equation (2.6)) is explicitly taken into account, can be written as

1 + ζ 2 ū(x) = S I U • n κ(r y ) -Bo h(r y ) dS y + (ζ -1) S I n • K • ū dS y + S W 1 U • T 1 • n dS y -ζ S W 2 U • T 2 • n dS y + V 1 U • F 1 dV y + ζ V 2 U • F 2 dV y (4.2)
where x locates a point on S I , n = n 12 and ζ = η 2 /η 1 is the viscosity ratio.

Here U and K are Green kernels for velocity and stress respectively and are given by [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF])

U (d) = 1 8π ( 1 d I + dd d 3 ), i = 1, 2 (4.3) K(d) = - 3 4π ( ddd d 5 ), i = 1, 2 (4.4)
where d = ȳx, ȳ(r y , z y ) is the integration point. The first term in the right hand side of equation (4.2) describes the flow contribution from interfacial tension and gravity, whereas the second term accounts for shear rates contrast on the interface. This term vanishes when ζ = 1. The third and fourth terms account for shear occurring on S W 1 and S W 2 as a result of the no-slip boundary condition. The two last terms account for the contribution of the flow due to the scattering force density. Unlike the other terms, they are evaluated on the volumes V 1 and V 2 respectively for fluid 1 and 2. Velocities on the interface as well as stress over all the boundaries S I , S W 1 and S W 2 are determined by solving the discrete form of this equation using a numerical procedure. This is performed by an analytical integration in the azimuthal direction [START_REF] Lee | The motion of a sphere in the presence of a deformable interface[END_REF][START_REF] Graziani | A boundary integral function method for axisymmetric viscous flows[END_REF]) reducing all surface boundaries S I , S W 1 , S W 2 into line boundaries Γ I , Γ W 1 , Γ W 2 and volumes V 1 and V 2 into surfaces S 1 and S 2 (see figure 2). The numerical procedure requires then the discretization of all line boundaries and surfaces with line segments and quadrilaterals respectively. Line discretization makes use of constant boundary elements, i.e. line segments with centered nodes. Elliptic line integrals resulting from the azimutal integration are evaluated using Gauss quadratures and power series expansions [START_REF] Bakr | The boundary integral equation method in axisymmetric stress analysis problems[END_REF]. Numerical integration on quadrilaterals is performed using two-dimensional Gauss quadratures [START_REF] Davis | Methods of numerical integration[END_REF]) and isoparametric transformations. The fluid-fluid interface is parameterized in terms of arc length and is approximated by local cubic splines, so that the curvature can be accurately computed. Distribution and number of points are adapted to the shape of the interface, so that the concentration of elements is higher in regions where the variation of curvature of the interface is larger.

The motion of the interface is followed using the kinematic condition (2.9) which is discretized using an explicit first-order Euler time scheme. A typical computation begins with a flat interface at rest. The laser beam is switched on at t=0, and the deformation of the interface begins. Computation stops when an equilibrium state is reached ( dh dt → 0). The time step is chosen to be about 50 times smaller than the reference time τ = L 1 /U σ .

Results

Single layer flow

In this section, we compare the steady velocity field in the bottom layer predicted using both the numerical simulation and the analytical one-layer model. For this comparison, we chose η 2 /η 1 = 10 -3 and ρ 2 << ρ 1 while f 2 = 0. Moreover, the influences of the layer aspect ratio R/L 1 and of the layer thickness to beam size ratio L 1 /ω 0 on the flow field and interface shape are discussed.

Flow pattern

Comparisons between the flow patterns predicted analytically and numerically are reported in figure 3. It shows the r and z dependences of the axial (figures The same remark holds for u r (figures 3(a,e-g)) at the horizontal mid-plane of the layer (z = -L 1 /2) and at the fluid-fluid interface (z = 0). A discrepancy can be identified at the bottom wall (z = -L 1 ) due to the no-slip condition (u r (z = -L 1 ) = 0) in the simulation while a free slip condition was chosen in the one-fluid model for simplicity purpose. For all radial velocities we can see that u r > 0 for z > -L 1 /2 and u r < 0 for z < -L 1 /2, which is consistent with the statement of a toroidal eddy centered vertically on the mid-plane of the fluid layer.

u z /U f -1 -0.8 -0.6 -0.4 -0.2 0 z/L 1 z/L 1 =-0.5 z/L 1 =-0.1 z/L 1 =-0.9 r/L 1 =0.2 r/L 1 =0.4 r/L 1 =0.8 (a) (b) (c) (d) (e) (f) (g) r/L 1 =0.1
Even though the one-layer model is a simplified version of the real situation, the comparison shows that this analytical model captures the essential trends of the flow. In the next section, we investigate the effects of the aspect ratio R/L 1 .

-

Influence of the container aspect ratio

Comparisons between the steady axial velocity profiles u z , predicted analytically and numerically, for different values of the aspect ratio R/L 1 , are reported in figure 4. It shows a good agreement between the simulations and the analytical model except for |r|/L 1 ≤ 0.1 i.e. within the beam |r| ≤ ω 0 , since L 1 /ω 0 = 10. This result was expected since the analytical solution shows a divergence of u z for r = 0 as a consequence of the line force distribution approximation. The left inset of figure 4 shows that u z scales like U f log(R/L 1 ) for narrow tubes (L 1 >> R) while it scales like U f in the thin films case (R >> L 1 ). The right inset of figure 4 shows that the radial extension of the flow r(u z = 0) scales like the smallest dimension of the container. In fact, when L 1 >> R, the flow is bounded by the lateral walls and consequently its radial extension scales like R. In the other limit, when R >> L 1 , the flow is unbounded in the r direction and mass conservation in the axisymmetric geometry produces toroidal eddies of circular section of radius L 1 . Deformation profiles predicted by the numerical simulations for different values of R/L 1 are presented in figure 5. The left inset shows that the tip deformation scales like L 1 when R >> L 1 . Similarly to the velocity field, the radial extension of the deformation scales like the smallest container dimension as shown by the right inset of figure 5. To conclude, in the limit of narrow beams (L 1 >> ω 0 ), the radial extensions of the flow and the deformation always scales like the smallest dimension of the layer. In the next section, the influence of the beam size on the flow and on the interface deformation is investigated.

Influence of the beam size

Comparisons between the steady axial velocity profiles u z , predicted analytically and numerically, for different values of the layer thickness to beam waist ratio L 1 /ω 0 , are reported in figure 6. It shall be noted that the analytical profile is unique since the solution does not depend on ω 0 . We observe that all numerical profiles coincide for |r| > ω 0 . In addition, we notice a disagreement for |r| < ω 0 , similarly to the results presented in figure 4. While the analytical prediction shows a divergence of u z as r → 0, the numerical simulation shows a saturation. This saturation always occurs at r ≃ ω 0 . We conclude, in the limit of thin films, that there are two characteristic dimensions for the axial velocity. When |r| < ω 0 , the axial velocity depends on ω 0 while its radial extension depends on the layer thickness L 1 (see right inset of figure 6). The left inset of figure 6 shows that u z scales like U f log(L 1 /ω 0 ) in the limit of thin films and narrow beams. A similar behaviour has been observed, in the Re << 1 limit, in a theoretical investigation of the acoustic streaming velocity field by Nyborg [START_REF] Hamilton | Nonlinear acoustics[END_REF]. Indeed, it was demonstrated that the maximum axial velocity of acoustic streaming increases logarithmically with the focal length (equivalent to L 1 ) of the acoustic beam propagating through an infinite fluid. Deformation profiles predicted by the numerical simulations for different values of L 1 /ω 0 are shown in figure 7. We conclude from the insets that both the deformation height and its radial extension scale like L 1 when L 1 >> ω 0 . The general conclusion of this section is that for thin films R >> L 1 and of narrow beams L 1 >> ω 0 , both the radial extension of the axial velocity and of the deformation scale like the layer thickness L 1 while their amplitudes scale respectively like U f and L 1 . These limits correspond to our experimental investigations reported previously [START_REF] Schroll | Liquid transport due to light scattering[END_REF][START_REF] Wunenburger | Fluid flows driven by light scattering[END_REF]. The following sections present the hydrodynamic interaction between the two layers with an investigation of the influence of the different aspect ratios R/L 1 , L 1 /L 2 and force and viscosity relative amplitude f 1 /f 2 and η 2 /η 1 .

Double-layer flow

Flow patterns

The flow patterns in wide containers (R > L 1 , L 1 = L 2 ) and long containers (R < L 1 , L 1 = L 2 ) are displayed in figure 8 and figure 9(a). We first observe in figure 8, nearly circular streamlines, which are the signature of toroidal flows also called eddies. These eddies exist at steady state because (i) the flow has to follow tangentially the interface due to the continuous axial forcing of the scattering force and (ii) mass conservation induces a circulation of the fluid. When R > L 1 , figure 8 clearly shows that the spatial extension of these eddies scales like L 1 as discussed in the previous sections (figure 4). When L 1 > R (figure 9(a)), the eddies invest the whole thickness of the layers while their width equals the container radius. To conclude, these results demonstrate that the axial dimension of the eddies is always L 1 while the radial one scales like the smallest dimension of the fluid layer.

Figures 9(b-c) display the flow pattern when L 1 and L 2 are different. We notice first that the maximum velocities are different in the layers despite an equal scattering force density. This is due to the dependence of the velocity with respect to the thickness of the fluid layer, which has been evidenced in the left inset of figure 6. It should be also noticed that in the fluid layer of aspect ratio L 1 /L 2 equal to 1, there is an eddy extending to the whole layer, while in the other layer of smaller height, a secondary eddy develops. In fact, because the eddy is larger in the fluid layer of largest thickness, it forces a motion in the other layer, as a result of momentum transfer by viscous shear through the interface, yielding to a secondary eddy. The scattering force exerted in both layers inducing two contrarotative eddies, any loss of symmetry of the flow domain induces the appearance of a secondary corotative eddy within one of the layers. Similar patterns can be obtained by imposing non equal scattering force densities in the two layers as illustrated in the next section.

Influence of the scattering force ratio

In the present section, we discuss the dependence of the flow and of the interface deformation on the magnitude of the scattering force. Figure 10 shows the variation of h max = h(r = 0) with the total scattering force f = f 1 + f 2 . We first notice that when the forces are equal in both layers f 1 = f 2 = f /2, the variation of h max is linear with respect to f . For [f 1 = 0, f 2 = f ] and [f 1 = f, f 2 = 0], we observe a small deviation from the linear behaviour, which is attributed to the dependence of h max to the layer thickness as shown in figure 7. Considering the case [f 1 = 0, f 2 = f ], the layer thickness above the hump (L 0 -h 2 ) is a little smaller than L 0 leading to a smaller deformation, while in the case [f 1 = f, f 2 = 0], the layer thickness below the hump (L 0 + h 1 ) is a little larger than L 0 leading to a larger deformation (see left inset of figure 7). The right inset of figure 10 shows that the variation of the difference of heights ∆h = h 1 -h 2 between the case [f 1 = f, f 2 = 0] and the case [f 1 = 0, f 2 = f ] is quadratic with respect to f . In the left inset of figure 10, we can see that the maximum velocities vary linearly with f , and that the same velocity is obtained in both configurations ([f 1 = f, f 2 = 0] and [f 1 = 0, f 2 = f ]). In the configuration f 1 = f 2 = f /2, the maximum velocity in each layer is half that obtained in the other configurations. Consequently, the effects of the flows on interface deformations are additive as each layer induces its own viscous stress on the interface. Considering the flow pattern, we observed in figure 8 that when Σ = f 1 /f 2 = 1, two co-rotating eddies are induced, both of them rotating clockwise, as a result of the direction of the scattering force density along the z axis. In this perfectly symmetric case, no secondary eddy is generated.

0 0.5 1 1.5 2 2.5 3 f/σ 0 0.05 0.1 0.15 h max /L 1 f 1 =f ; f 2 =0 f 2 =f ; f 1 =0 f 1 =f 2 =f/2 0 1 2 3 4 f/
In figures 11, 12 and 13, corresponding to Σ = 1, a secondary eddy rotating counter-clockwise develops in the fluid layer where the amplitude of the scattering force is the smallest. This secondary eddy, whose extension increases with Σ, is contra-rotative compared to the main eddies. Finally, in figure 14, corresponding to [f 1 = f, f 2 = 0], we observe two contrarotating eddies, the one in the bottom layer being determined by the direction of propagation of the laser beam while the one in the top layer, which is of weaker intensity, is exclusively induced by viscous shear through the interface as a result of zero scattering force density, f 2 = 0.

To summarize, when f σ, there is no difference whether the scattering force is applied in the top or in the bottom layer. When f > σ, a small difference between the deformations can be noticed and is attributed to the resulting difference of thickness of the scattering layer. The flow patterns showed that a difference of scattering force between the layers leads to the emergence of a secondary eddy in the layer where the force is the smallest.

Up to now, we have investigated flow patterns for ζ = η 2 /η 1 = 1 in order to analyze the influence of the scattering force amplitude and of the different aspect ratios of the container. However, since viscosity affects flow velocities, the question remains on whether it modifies the two-layer flow pattern driven by light scattering. This is the purpose of the following section.

Influence of the viscosity ratio

In figure 15, the variations of the maximum velocity

u max i in fluid layer i is represented versus ζ = η 2 /η 1 while keeping Σ = f 1 /f 2 = 1 and L 1 = L 2 .
As we set η 1 to a constant value, there is no variation of the velocity magnitude in layer 1, while the velocity in layer 2 scales like 1/ζ = η 1 /η 2 . This demonstrates that the dimensional scaling law established for a single layer flow

u i ∝ F i ω 2 0 /η i (f i = F i πω 2 0
) is valid in each fluid layer as ω 0 is a characteristic length for the axial velocity near the z-axis (see comments of figure 6). Moreover, the inset of figure 15 shows that h max does not depend on the viscosity ratio. This can be explained by recalling that T zzi is proportional to η i u i while u i is proportional to 1/η i making T zzi and hence h max independent of the fluid viscosities. In order to get a general understanding of the location of the secondary eddy, a parametric study of the flow pattern was performed by varying Σ and ζ in the set (Σ, ζ) ∈ {0.01; 1; 100} 2 while keeping the total scattering force f = f 1 + f 2 con- stant. Table 1 summarises these results, indicating whether the secondary eddy is in the bottom layer (1), top layer (2) or nonexistent (0). Table 1 shows that the location of the secondary eddy is independent of ζ and is always located in the layer were the scattering force density is the smallest. This can be explained by noting that, since F 1 ∝ η 1 u 1 /ω 2 0 and F 2 ∝ η 2 u 2 /ω 2 0 , one has u 1 /u 2 ∝ Σζ. This scaling law is confirmed by our computational results, as illustrated by the two following cases. In the first case, Σ = 10 and ζ = 0.1, we found that the maximum velocities in layers 1 and 2 are u max 1 /U σ ≃ u max 2 /U σ ≃ 0.012. In the second case, Σ = 0.1 and ζ = 0.1, we found u max 1 /U σ ≃ 0.00114 and u max 2 /U σ ≃ 0.114, in agreement with the prediction u 1 /u 2 = 0.01. The viscosity ratio only determines the velocity magnitude in each fluid layer. This observation is confirmed by a simple model of parallel flow which is presented in the following section.

Lubrication model

The location of the secondary eddy is predicted by making use of a lubrication model assuming R >> L i >> ω 0 and derived from a two-dimensional cartesian formulation of the Stokes equation. For the sake of simplicity, we consider a twodimensional flow of velocity u = (u x , 0, u z ) in cartesian coordinates (x, y, z). Since we are only interested in the distribution of eddies i.e. in the z dependence of u x , we consider layers of infinite horizontal extent and a parallel flow u = u x (z)e x . Therefore far enough from the z-axis along which the forcing is exerted (x >> 0), we can write the Stokes equation in both layers, as follows

(a) 0 = - ∂p i ∂x + η i ∂ 2 u ix ∂z 2 ; (b) 0 = - ∂p i ∂z ; i = 1, 2 (5.1)
We can deduce from (5.1 (b)) that the flow is driven by a horizontal pressure gradient.

Differentiating equation ( 5.1 (a)) with respect to z, we find that

0 = ∂ 3 u ix ∂z 3 ; i = 1, 2 (5.2)
This leads to the following expression of the velocity in each layer

u ix (z) = a i z 2 + b i z + c i ; i = 1, 2 (5.3)
where a i , b i , c i , i = 1, 2 are six numerical constants to be determined by using the six following conditions

u 1x (z = 0) = u 2x (z = 0) (5.4) η 1 du 1x dz (z = 0) = η 2 du 2x dz (z = 0) (5.5) 0 -L 1 u 1x (z)dz = 0 (5.6) L2 0 u 2x (z)dz = 0 (5.7) dp 1 dx = -F 1 (5.8) dp 2 dx = F 2 (5.9)
The first two conditions (5.4) and (5.5) simply express the continuity of velocities and shear stress along the interface. Conditions (5.6) and (5.7) express mass conservation in each layer. Conditions (5.8) and (5.9) express the fact that the pressure gradient driving the flow is related to the scattering force driving the motion near the axis as explained hereafter. Because the eddies always scale like L i when R >> L i , we consider that ∂p i ∂x ∝ ∂p i ∂z ∝ p i /L i . In addition, the Stokes equation near the z-axis leads to ∂p i ∂z ∝ η i ∂ 2 u iz ∂x 2 ∝ F i and finally to ∂p i ∂x ∝ p i /L i ∝ F i . One can also notice a (-) sign on equation (5.8). This is justified by the fact that when the system is symmetric (i.e. F 1 =F 2 , η 1 = η 2 , L 1 = L 2 ), we assume that there is no motion of the interface along the x axis, which leads to opposite pressure gradients and velocities near the interface (u 1x (z = 0 -) > 0 and u 2x (z = 0 + ) < 0). The solution approximated the flow without satisfying the no-slip boundary conditions at the top and bottom walls. The conditions imposed to the two-fluid, parallel flow lead to a 1 = -F 1 2η 1 (5.12)

c 1 = c 2 = F 1 6η 1 L 2 1 - L 1 ( F 1 L 2 1 3η 1 + F 2 L 2 2 3η 2 ) 2(L 1 + L 2 /ζ) (5.10) b 1 = ζb 2 = - F 1 L 2 1 3η 1 + F 2 L 2 2 3η 2 L 1 + L 2 /ζ (
a 2 = F 2 2η 2 (5.13)
Figure 16 illustrates four interesting cases highlighting the emergence or not of a secondary eddy in one of the layers. In figure 16(a), the system is symmetric and u x (z) changes sign only once in each layer. This means that only one eddy exists in each layer. In case (b), the magnitude of the scattering force density in the top layer is half that in the bottom one. In the top layer, u x (z) changes sign twice (at z = 0.79 and z = 0.21). This means that two contra-rotative eddies exist in the top layer. In case (c), the fluid in the bottom layer is 5 times more viscous than the one in the top layer. We can notice that although the velocities in the bottom layer are much smaller, there is no secondary eddy. Finally, in case (d) the bottom layer is 1.5 times thicker than the top layer. A secondary flow appears in the layer of smallest thickness, since u x (z) changes sign twice at z = 0.73 and z = 0.14. As a conclusion, even though very simple, this parallel flow model gives a qualitative prediction of the location of the secondary eddy confirming the numerical simulations of the axisymmetric flow.

General discussion

Optical streaming represents a new way to drive fluid flows by light. It can be useful in situations where the direction and magnitude of fluid flows need to be controlled. The only condition for inducing these optical streaming flows is to use non absorbing turbid fluids. Any viscous fluid can be considered, while the strength of the optical streaming will strongly depend on the size and nature of the scatterers which could be solid or liquid nanosuspensions (see [START_REF] Wunenburger | Fluid flows driven by light scattering[END_REF]; [START_REF] Schroll | Liquid transport due to light scattering[END_REF] for an example of density fluctuations in nearcritical binary mixtures acting as liquid scatterers) and generally speaking any fluid presenting refractive index inhomogeneities at a length scale much smaller than the optical wavelength. Microfluidics is of particular interest for application of optical streaming since the small dimensions involved are compatible with the laser focusing. If a fluid-fluid interface is present in these experiments (either microfluidics coflow experiments or digital microfluidics), we expect the emergence of secondary eddies depending on the thicknesses of the fluid layers and on the scatterers present in these fluids. In addition, if the interfacial tension is sufficiently low, this could also lead to interface deformation. While a viscosity contrast will not influence the location of secondary flows, a significant impact of this contrast on flow magnitude is expected since velocity is inversely proportional to viscosity. Many applications of optical streaming can thus be devised. For instance, flow control at very small flow rates is always difficult in microchannels as syringe pumps operate close to their mechanical limits; flow rate changes are also always accompanied by complex transients. In these regimes, the fine and smooth control of viscous flows by optically tunable scattering forces is probably an interesting alternative for turbid fluids. Micromixing is also very challenging in classical microchannels [START_REF] Nguyen | Micromixersa review[END_REF] because flows are mostly laminar and thus constitute a major drawback for efficient chemical reactions [START_REF] Baroud | Reaction-diffusion dynamics: Confrontation between theory and experiment in a microfluidic reactor[END_REF]. One route for active mixing consists in superimposing pressure perturbations to the main stream [START_REF] Lee | Chaotic mixing in electrokinetically and pressure driven micro flows[END_REF]. Scattering forces could replace pressure actuation if fluids are turbid or seeded with appropriate scatterers.

Conclusion

In this work, we investigated flow production by light scattering in a two-fluid system involving turbid fluids. We compared the flow pattern predicted numerically to the analytical one-fluid model prediction. This simplified analytical description predicts most of the characteristic features of both the induced flow and the resulting interface deformation. Due to the large scale of induced flows, we investigated the geometrical effects of the container. We showed that the velocity field and the interface deformation scale like the smallest dimension of the container. Moreover a variety of flow patterns were illustrated when breaking the system symmetries. For instance, in absence of scattering force intensity contrast, secondary eddies always emerge in the thinnest layer while in absence of layer thickness contrast, secondary eddies emerge in the fluid layer where the scattering force density is the weakest. On the contrary the variation of the viscosity ratio only affects velocities within the fluids while neither the interface deformation nor the flow pattern at steady state are influenced. Finally, a lubrication model was used to explain the emergence and to predict the location of these secondary flows in the layers. Generally speaking, light-induced flows exist whenever isotropic scattering fluids are at work; fluids near a second-order phase transition can be an interesting example to perform experiments. Consequently, we have proposed and wellidentified a general non dissipative way to produce bulk flow with light that goes one step forward in the exploration of the coupling between light and hydrodynamics, so called "optohydrodynamics".

8. Appendix : Comparison with experiments : Interface deformation in the presence of radiation pressure In this appendix, we show comparisons between numerical, analytical and experimental steady profiles of interface deformations. In experiments involving turbid two-phase systems like the one used in [START_REF] Schroll | Liquid transport due to light scattering[END_REF], the interface is made visible thanks to a refractive index contrast. In addition to the hydrodynamics stress due to the optical streaming, radiation pressure is exerted on the interface (see [START_REF] Chraibi | Optohydrodynamics of soft fluid interfaces: Optical and viscous nonlinear effects[END_REF]). In this appendix, both stresses are taken into account in the numerical simulation. Note that in [START_REF] Schroll | Liquid transport due to light scattering[END_REF], only optical streaming effects were taken into account while in [START_REF] Chraibi | Optohydrodynamics of soft fluid interfaces: Optical and viscous nonlinear effects[END_REF], only radiation pressure effects were modeled. In order to model radiation pressure effects, we have to add the following term (which is presented in details in [START_REF] Chraibi | Optohydrodynamics of soft fluid interfaces: Optical and viscous nonlinear effects[END_REF]) in the right side of equation (2.6) Π(r) = I(r) c n 2 cos θ i (2(1 -δ) cos θ i -Υ(θ i , θ t )((1 -δ) cos θ i + cos θ t )) (8.1)

The radiation pressure is due to the jump in photon linear momentum (proportional to the refractive index n i ) when propagating from fluid 1 to fluid 2. Υ(θ i , θ t ) is the Fresnel transmission coefficients of energy fluxes for circularly polarized beams (Chraibi et al. 2008a)

Υ(θ i , θ t ) = 2(1 -δ) cos θ i cos θ t ((1 -δ) cos θ i + cos θ t ) 2 + 2(1 -δ) cos θ i cos θ t (cos θ i + (1 -δ) cos θ t ) 2 (8.2)
where θ i and θ t are the incidence and transmission angles at the fluid interface respectively. They are such that θ i = arctan( dh dr ) and θ t = arcsin((1 -δ) sin θ i ) where δ = (n 2 -n 1 )/n 2 is the relative refractive index contrast. The fluid used in this experiment is a near-critical binary mixture in the twophase regime that undergoes refractive index fluctuations of characteristic size ξ -, which increases as the critical point is neared. When the difference between the sample temperature T and the critical temperature T c is large, ξ -vanishes and each phase at coexistence can be considered as almost homogeneous [START_REF] Casner | Giant deformations of a liquid-liquid interface induced by the optical radiation pressure[END_REF][START_REF] Casner | Asymmetric optical radiation pressure effects on liquid interfaces under intense illumination[END_REF] for the light-matter interaction. However, when T is close to T c , ξ -increases up to tens of nanometers at T-Tc=0.1K and the sample significantly scatters light, producing strong linear momentum transfer to the fluid. Details on the experiments will not be presented here, however one can refer to [START_REF] Schroll | Liquid transport due to light scattering[END_REF], [START_REF] Wunenburger | Fluid flows driven by light scattering[END_REF] and [START_REF] Chraibi | Optohydrodynamics of soft fluid interfaces: Optical and viscous nonlinear effects[END_REF] for more details. The profiles are analytically predicted using the stream function formulation developed in section 3 all along with the numerical solution of equation (3.16) (considering h ′ (r = 2ω 0 ) = 0 as observed in the experiments). Figure 17 presents comparisons between two-fluids numerical results, one-fluid analytical predictions and experimental steady interface profiles for two different experimental conditions (T-T c =1 K, ω 0 =5.3 µm, and T-T c =1.5 K, ω 0 =7.5 µm) and at four different values of the laser beam power. A very good agreement is obtained between the experimental data and the predictions at each power value P = 154, 308, 462 and 616 mW. Insets of figure 17 illustrate the linear power dependence of the deformation height due to the steady scattering flows disturbed by opposite radiation pressure effects. All these results validate both the numerical resolution and the theoretical model, which both represent important tools to describe these original physical phenomena.

Figure 2 .

 2 Figure 2. Sketch of the flow domain composed of the two fluids and the liquid interface initially at rest.

  r, z = -L 1 /2) = 0, the corrected pressure field can be deduced from equation (3.2) by noticing that ∂p ∂r =

Figure 3 .

 3 Figure3. Comparison between numerical (symbols) and analytical (lines) velocity fields in the bottom layer. The left figures show the variation of the reduced velocity components ur/U f and uz/U f with respect to r/L1 (U f = f1/η1). The right ones present the variation of ur/U f and uz/U f with respect to z/L1. f1/σ = 0.157, f2 = 0, η2/η1 = 10 -3 , R/L1 = 5, L1 = L2 and L1/ω0 = 10. In the numerical resolution Bo = 5 and therefore h(r = 0)/L1 = 0.008.

Figure 4 .

 4 Figure 4. Axial velocity for several values of R/L1. f1/σ = 0.157, f2 = 0, η2/η1 = 10 -3 , L1 = L2 and L1/ω0 = 10. In the numerical simulation Bo = 5. Open symbols represent the numerical simulation while continuous lines represent the analytical resolution. Left inset shows the variation of uz(r = 0.05)/U f versus R/L1. Right inset shows the variation of r(uz = 0)/L1 versus R/L1.

Figure 5 .

 5 Figure 5. Numerical deformation profiles for several values of R/L1. f1/σ = 0.157, f2 = 0, η2/η1 = 10 -3 , Bo = 5, L1 = L2 and L1/ω0 = 10. Left inset shows the variation of h(r = 0)/L1 versus R/L1. Right inset shows the variation of r(h = 0)/L1 versus R/L1.

Figure 6 .

 6 Figure 6. Axial velocity for several values of L1/ω0. f1/σ = 1.57, f2 = 0, η2/η1 = 10 -3 , L1 = L2 and R/L1 = 5. In the numerical simulation Bo = 5. Open symbols represent the numerical simulation while the continuous line represents the analytical resolution. Left inset shows the variation of uz(r = 0.02)/U f versus L1/ω0. Right inset shows the variation of r(uz = 0)/L1 versus L1/ω0. Beam waist location is shown by the dashed lines for each profile.

Figure 7 .

 7 Figure 7. Numerical deformation profiles for several values of L1/ω0. f1/σ = 1.57, f2 = 0, η2/η1 = 10 -3 , Bo = 5, L1 = L2 and R/L1 = 5. Left inset shows the variation of the tip height h(r = 0)/L1 versus L1/ω0. Right inset shows the variation of the radial extension of the deformation r(h = 0)/L1 versus L1/ω0.

Figure 8 .

 8 Figure 8. Steady flow pattern for R/L1 = 5 and L1/ω0 = L2/ω0 = 30, Σ = f1/f2 = 1, ζ = η2/η1 = 1, Bo = 10 and f1/σ = 0.3. Velocity magnitude is reduced by Uσ = σ/η1.

Figure 9 .

 9 Figure 9. Steady flow patterns. (a) R/L1 = 2/3 and L1/ω0 = L2/ω0 = 30 (b) R/L1 = 3, L1/ω0 = 10 and L2/ω0 = 30 (c) R/L1 = 1, L1/ω0 = 30 and L2/ω0 = 10. In all cases f1/σ = 0.3, Σ = f1/f2 = 1, ζ = η2/η1 = 1 and Bo = 10. Velocity magnitude is reduced by Uσ.

Figure 10 .

 10 Figure 10. Hump height variation against the total scattering force f when [f1 = f, f2 = 0], [f1 = 0, f2 = f ] and f1 = f2 = f /2 (R/L1 = 3.33, L1/ω0 = L2/ω0 = 30, Bo = 45 and η1 = η2). The dashed line is a linear regression for f1 = f2 = f /2. The left inset presents the maximum fluid velocity in each configuration while the right inset shows the variation of the reduced difference of heights ∆hmax/L1 between the cases [f1 = f, f2 = 0] and [f1 = 0, f2 = f ] against f /σ.

Figure 11 .

 11 Figure 11. Steady flow pattern for R/L1 = 5 and L1/ω0 = L2/ω0 = 30, Σ = f1/f2 = 2, ζ = η2/η1 = 1, Bo = 10 and f1/σ = 0.3. Velocity magnitude is reduced by Uσ.

Figure 12 .

 12 Figure 12. Steady flow pattern for R/L1 = 5 and L1/ω0 = L2/ω0 = 30, Σ = f1/f2 = 4, ζ = η2/η1 = 1, Bo = 10 and f1/σ = 0.3. Velocity magnitude is reduced by Uσ.

Figure 13 .

 13 Figure 13. Steady flow pattern for R/L1 = 5 and L1/ω0 = L2/ω0 = 30, Σ = f1/f2 = 1/4, ζ = η2/η1 = 1, Bo = 10 and f1/σ = 0.3. Velocity magnitude is reduced by Uσ.

Figure 14 .

 14 Figure 14. Steady flow pattern for R/L1 = 5 and L1/ω0 = L2/ω0 = 30, f1/σ = 0.3, f2 = 0, ζ = η2/η1 = 1, Bo = 10. Velocity magnitude is reduced by Uσ.

Figure 15 .

 15 Figure 15. Variation of the maximum velocities in each fluid layer as a function of the viscosity ratio ζ = η2/η1. η2 is varied whereas η1 is kept constant (Σ = f1/f2 = 1, f1/σ = 0.3, R/L1 = 3.33, L1/ω0 = L2/ω0 = 30, Bo = 45). The inset shows the variation of hmax/L1 with ζ.

Figure 16 .

 16 Figure 16. Velocity ux(z) predicted by the cartesian model. (a) Symmetric case: Σ = 1, ζ = 1 and L1/L2 = 1. (b) Contrast of magnitude of scattering force densities: Σ = 2, ζ = 1 and L1/L2 = 1. (c) Contrast of viscosities: Σ = 1, ζ = 0.2 and L1/L2 = 1. (d ) Contrast of layer thicknesses: Σ = 1, ζ = 1, L1/L2 = 1.5.

Table 1 .

 1 Secondary eddy location for (Σ, ζ) ∈ {0.1; 1; 10} 2 . Bottom layer, top layer or inexistent secondary eddy are respectively marked as 1, 2 and 0.

† E-mail address: h.chraibi@loma.u-bordeaux1.fr


Igal Milchtaich

Representation of Finite Games as Network Congestion Games

Keywords: network games, congestion games, potential games I

Weighted network congestion games are used for modeling interactions involving finitely many non-identical users of network resources, such as road segments or communication links. In spite of their special form, these games are not fundamentally special: every finite game can be represented as a weighted network congestion game. The same is true for the class of unweighted network congestion games with player-specific costs, in which the players differ in their cost functions rather than their weights. The intersection of the two classes consists of the unweighted network congestion games. These games are special: a finite game can be represented in this form if and only if it is an exact potential game.

INTRODUCTION

In a (finite) congestion game, finitely many players share a finite set of resources but may differ in which resources they are allowed to use. Specifically, each player's set of strategies is a particular collection of nonempty subsets of . The player's payoff from using a strategy is the negative of the total cost of using the resources included in the strategy. The cost of a resource depends only on its identity and on the number of users. The cost does not necessarily increase with congestion, and it may be negative, which means that using the resource contributes positively to the payoff. Rosenthal [START_REF] Rosenthal | A class of games possessing pure-strategy Nash equilibria[END_REF] showed that every congestion game admits an exact potential, which is a function on strategy profiles that exactly mirrors the players' incentives to unilaterally change their strategies. Whenever a single player moves to a different strategy, his gain or loss is equal to the corresponding change in . Monderer and Shapley [START_REF] Monderer | Potential games[END_REF] proved the converse: essentially, only congestion games are exact potential games. More precisely, every finite game that admits an exact potential can be represented as a congestion game. One of this paper's findings strengths Monderer and Shapley's result by showing that an exact potential game can always be represented as a particular kind of congestion game, namely, an unweighted network congestion game.

A restriction or expansion of the meaning of 'congestion game' potentially has the same effect on the class of representable games. Examples of restriction are: congestion games with nondecreasing cost functions, in which increasing congestion never makes users better off; singleton congestion games, in which each strategy includes a single resource; and network congestion games, in which resources are represented by edges in a graph and strategies correspond to routes, which are paths in the graph that connect the player's designated origin and destination vertices. (The second restriction is a special case of the third. It corresponds to a parallel network, which is one with only two vertices.) Examples of extensions are: congestion games with player-specific costs (or payoffs [START_REF] Milchtaich | Congestion games with player-specific payoff functions[END_REF]), where players are differently affected by congestion; and weighted congestion games, in which their contributions to it (the players' "congestion impacts") differ.

This paper shows that both weighted network congestion games and unweighted network congestion games with playerspecific costs are actually capable of representing all finite games. Both representations use only nondecreasing cost functions and two-terminal, or single-commodity, networks, in which all players' routes start and terminate at the same origin and destination vertices. Thus, although the definitions of the two kinds of network congestion games involve quite specific structures, the games themselves are not in any way special. For unweighted network congestion games, which are simultaneously weighted network congestion games and unweighted network congestion games with player-specific costs, this is not so. As indicated, these games are special in that they represent (all) exact potential games.

The potential significance of these findings lies in the information that the representation as a network congestion game provides about the represented game. An example of such information is existence of a pure-strategy equilibrium [START_REF] Milchtaich | Network topology and equilibrium existence in weighted network congestion games[END_REF]. For certain network topologies, the existence of equilibrium is guaranteed regardless of the exact functional form of the cost functions, the number of players and (if applicable) their weights. Thus, any finite game that can be represented as a network congestion game on such a network has at least one pure-strategy equilibrium.

II. PRELIMINARIES

A. Game Theory

A finite (noncooperative) game has a finite number of players whose strategy sets are finite. A strategy profile in assigns a strategy to each player . Two games and with identical sets of players are isomorphic [START_REF] Monderer | Potential games[END_REF] if for each player there is a one-to-one correspondence between 's strategy sets in and , such that each strategy profile in yields the same payoffs to all players as the corresponding strategy profile in .

Essentially, isomorphic games are just alternative representations of a single game. Two games and with identical sets of players and respective strategy sets are similar if for each player the difference between the payoffs in and does not depend on the player's own strategy, and can thus be expressed as a function of the other players' strategies.

A game is an exact potential game [START_REF] Monderer | Potential games[END_REF] if it is similar to some game in which all the players have the same payoff function , which is said to be an exact potential for .

B. Graph Theory

An undirected multigraph consists of a finite set of vertices and a finite set of edges. Each edge joins two distinct vertices, and , which are referred to as the end vertices of . Thus, loops are not allowed but more than one edge can join two vertices. An edge and a vertex are incident with each other if is an end vertex of . A (simple) path of length is an alternating sequence of vertices and edges , beginning and ending with vertices, in which each edge is incident with the two vertices immediately preceding and following it and all the vertices (and necessarily all the edges) are distinct. Every path traverses each of its edges in a particular direction: from the end vertex that precedes in the path to the vertex that follows it.

A two-terminal network, or simply network, is an undirected multigraph together with a distinguished ordered pair of (distinct) terminal vertices, the origin and the destination , such that each vertex and each edge belongs to at least one path in which the first vertex is and the last vertex is . Any path with these first and last vertices is referred to as a route in .

Two networks and may be connected in parallel if they have the same origin and the same destination but no other common vertices or edges, and in series if they have only one common vertex which is the destination in and the origin in . In both cases, the set of vertices and the set of edges in the resulting network are the unions of the corresponding sets in and , and the origin and destination are those in and , respectively (as well as in and , respectively, in the case of connection in parallel). The connection of an arbitrary number of networks in parallel or in series is defined recursively.

C. Network Congestion Games

A weighted network congestion game on a (two-terminal 1 ) network is a finite, -player game that is defined as follows. First, an allowable direction, a (possibly, empty) set of allowable users and a nondecreasing cost function are specified for each edge in , such that each edge is traversed in the allowed direction by at least one route and each player has at least one allowable route, that is, a route in that includes only edges that is allowed to use and traverses them in the allowed direction. The strategy set of 1

The assumption of a single origin-destination pair may be viewed as a normalization. Any weighted network congestion game on a multicommodity network, which has multiple origin-destination pairs, may also be viewed as a game with a single such pair. In that game, each terminal vertex is incident with a single allowable edge for each player, which joins it with the player's corresponding terminal vertex in the original game.

each player is the collection of his allowable routes. Second, a weight is specified for each player , which represents the player's congestion impact. 2 The total weight of the players whose chosen route includes an edge is the flow (or load) on . The (not necessarily negative) cost of for each of its users is . A player's payoff in the game is the negative of the total cost of the edges in his route.

A weighted network congestion game is called an unweighted network congestion game if the players' weights are all identical and equal to . The equality entails, in particular, that the cost of an edge is not affected by the identities of its users but only by their number. A generalization that allows for a dependence of the cost on the user's identity is unweighted network congestion game with player-specific costs. In such a game, each edge is associated with a (player-specific) nondecreasing cost function for each player , and its cost for that player is

, where (the flow) is the total number of players using .

III. PRESENTATION RESULTS

Every unweighted network congestion game is in particular a congestion game in the sense of belonging to the class of games presented by Rosenthal [START_REF] Rosenthal | A class of games possessing pure-strategy Nash equilibria[END_REF] and studied by Monderer and Shapley [START_REF] Monderer | Potential games[END_REF] (see the Introduction). In fact, as the following theorems show, the class of unweighted network congestion games essentially coincides with that of all congestion games. Weighted network congestion games and unweighted network congestion games with player-specific costs are generally not congestion games (in the above sense). The next theorem shows that both classes essentially coincide with that of all finite games. Figure 1.

Representations of a exact potential game (a) as a weighted network congestion game (b), with weights and , and as an unweighted network congestion game with playerspecific costs (c). Dotted, dashed and solid edges are allowable to player 1, player 2 and both players, respectively. The allowable directions are indicated where required. All relevant costs other than those specified are zero. A player's payoff is the negative of his total cost.

2

In certain contexts, it may be desirable to require the following (weak) connection between the players' weights and the cardinality of their strategy sets: For all and with , . In this paper, this requirement is not needed, but adding it would not affect any of the paper's results.

Theorem 1. Every finite game is isomorphic both to a weighted network congestion game and to an unweighted network congestion game with player-specific costs . is isomorphic to an unweighted network congestion game 3 if and only if it is an exact potential game.

By the first part of the theorem, every finite game can be represented as a network congestion game where the players differ in their weights and as one where they differ in their costs functions. The existence of a representation of the latter kind was first indicated by Monderer [START_REF] Monderer | Multipotential games[END_REF]. However, Monderer's definition of network congestion game is different, and significantly less restrictive, than in this paper. The second part of Theorem 1 strengths a well-known result of Monderer and Shapley [START_REF] Monderer | Potential games[END_REF]. Their result is the equivalence of (ii) and (iii) in the following extension of that part of the theorem.

2.

For every finite game , the following conditions are equivalent:

(i)
is isomorphic to an unweighted network congestion game.

(ii) is isomorphic to a congestion game.

(iii) is an exact potential game.

A finite game obviously has more than a single pair of representations as in Theorem 1. The "canonical" games and constructed in the proof of the theorem, which share the same network , are just one such pair. Other representations may be preferable in that they use a simpler network. An example is shown in Figure 1(c). That representation of the exact potential game in Figure 1(a) as an unweighted network congestion game with player-specific costs uses the Wheatstone network, which has fewer edges than the network constructed in the proof of Theorem 1 (which, for all games, is the one in (b)). Both networks are considerably simpler than the one constructed in the proof of Theorem 2. For an example of a simple representation of a specific variety of congestion games as weighted network congestion games, see [START_REF] Milchtaich | Weighted congestion games with separable preferences[END_REF].

Proof of Theorem 1. The following proof concerns the first part of the theorem. The second part is dealt with as part of the proof of Theorem 2 below.

Suppose that the number of players in and the cardinality of the largest strategy set are both at least two (otherwise the assertion is trivial). Without loss of generality, it may be assumed that, for all , player 's number of strategies is at least as great as that of (otherwise take 'player 1' below to mean the player with the largest number of strategies, 'player 2' the player with the second-largest number, and so on). It is desirable to have the stronger properly that these numbers are actually equal: all players have strategies. Such equality may be achieved by (temporarily) increasing the number of strategies of one or more players . To this end, one of 's strategies needs to be replicated, such that choosing the original strategy or any of its replicas has the same effect on 's payoff and those of the other players. Each player's strategies can now be indexed from to , thus identifying the 3 This condition can be expressed as the requirement that .

collection of all strategy profiles with the product set . Order this set in the following way:

where the order of the elements represented by the second ellipsis mark is immaterial. With each strategy profile associate two vertices and and an edge that joins them. The edge is directed from to and is public, that is, allowable to all players. Next, for each player and , consider all strategy profile with and list them according to their order in [START_REF] Milchtaich | Congestion games with player-specific payoff functions[END_REF]. For each pair and of successive entries in this list, add an edge that joins and , is directed from to , and is private, that is, allowable to player only. Finally, identify all vertices of the form , where is one of the first elements in [START_REF] Milchtaich | Congestion games with player-specific payoff functions[END_REF], and denote this single vertex by . Do the same for all vertices of the form , where is one of the last elements in [START_REF] Milchtaich | Congestion games with player-specific payoff functions[END_REF], and denote the result by . These terminal vertices, together with the other vertices and edges specified above, constitute a network , with specified direction and set of allowable users for each edge. (For games, this is the network depicted in Figure 1(b).) For each player , each allowable route in corresponds to some strategy in . Specifically, includes all edges with , alternating with private edges. Different allowable routes to player have no shared edges, and no shared vertices other than the terminal ones.

The unweighted network congestion game with playerspecific costs is defined by assigning the following cost functions to each player . For a public edge , corresponding to a strategy profile , where is player 's payoff function in and is any number equal to or greater than that function's maximum. For a private edge allowable only to player , . As explained above, strategy profiles in are in a one-to-one correspondence with allowable route choices in . The routes that correspond to a strategy profile are such that exactly one edge, namely , is used by all players. Therefore, for each player , only and private edges make a nonzero contribution to the cost. By (2), the total cost is . Hence, player 's payoff is , as in .

The weighted network congestion game is defined by, first, attaching the weight to each player ( ). Thus, the weight uniquely identifies the player, and is less than the total weight of any two players, which is or greater. Second, the cost functions are defined as follows. For a public edge , corresponding to a strategy profile ,

where is player 's payoff function in and is any number large enough to make the cost functions of all public edges nondecreasing. In this definition, is the collection of all partial strategy profiles , is defined as if the indicated equality holds and otherwise, and . For a private edge , which is allowable only to a single player , the cost is given by These "private" costs cancel the total contribution of the second term in the definition of the "public" costs, which appears only in the public edges in which player is the sole user. The total contribution of the corresponding first term, which depends on the strategy profile , can be computed as follows (to enhance readability, the computation is shown for ):

Again, the payoff is , as in .

To complete the proof of the first part of the theorem it only remains to dispose of the spurious strategies (potentially) introduced by replication at its beginning. These strategies may be eliminated simply by deleting or disallowing the use of all the private edges that belong to the corresponding routes.

Proof of Theorem 2. The fact that every congestion game, and in particular every unweighted network congestion game, is an exact potential game is well known [START_REF] Monderer | Potential games[END_REF], [START_REF] Rosenthal | A class of games possessing pure-strategy Nash equilibria[END_REF]. The reverse implications are proved below, using a variant of the proof in [START_REF] Monderer | Potential games[END_REF].

Suppose that has an exact potential . As in the proof of Theorem 1, and for similar reasons, it is sufficient to consider the case in which the number of players is at least and the strategy set of each player includes at least two elements. Strategy profiles in are naturally identifiable with subsets of the disjoint union (where 'disjoint' means that strategies of different players are viewed as distinct elements): is identified with the set . Similarly, a partial strategy profile , which is obtained from a strategy profile by ignoring the coordinate corresponding to a particular player , is identified below with the set . Player (who, because of the assumption of at least two strategies for each player, is the only player with

) is said to own . A congestion game isomorphic to is defined as follows. Each strategy profile in and each partial strategy profile are viewed as resources. (Thus, a resource is identified with a specific subset of .) For each player and strategy of that player in , the corresponding strategy in the congestion game is the selection of all the resources that include . These resources are of three kinds: all partial strategy profiles owned by , all partial strategy profiles owned by another player but containing , and all strategy profiles containing . The resources of the first kind are included also in every other strategy of player in the congestion game but those of the other two kinds are exclusive. It follows that player is the only user of (the resource identified with) a partial strategy profile if and only if (a) he owns it, i.e., it is of the form , and (b) the strategy of every other player is not (that corresponding to) , player 's coordinate in . The cost for a player of using a resource (which is the negative of the resource's contribution to the player's payoff) is defined as follows. For (the resource identified with) a strategy profile , the cost is if all players use the resource and otherwise. For a partial strategy profile , the cost is if no one else uses the resource (which, as indicated above, is possible only if the player considered is himself) and otherwise. In (3), is player 's payoff function in , is the collection of all strategy profiles, and is defined as if the indicated equality holds and otherwise. Note that a sufficient condition for the cost functions to be nondecreasing is that both and (3) are negative. Without loss of generality, this condition may be assumed to hold. It is not difficult to show that subtracting any constant from leaves it an exact potential function, and subtracts a certain positive multiple of from (3).

To prove that the congestion game defined above is isomorphic to it has to be shown that, for every player and strategy profile , the total cost of 's route is , the negative of his payoff in . Two kinds of resources make a nonzero contribution to the cost: the single resource corresponding to the strategy profile , which contributes , and those corresponding to partial strategy profiles owned by in which the strategy of every player is not . Assuming, for readability, that , the total contribution of the latter is The last equality uses the fact that, by definition of exact potential, for every player and strategy profile the difference does not depend on 's strategy . This completes the proof of the isomorphism between and the (3) congestion game, which corresponds to the implication (iii) (ii).

The above proof also constitutes the main step in the proof of the stronger assertion (iii) (i). To complete the proof, the congestion game defined above needs to be turned into an unweighted network congestion game. To this end, the resources in the congestion game need to be viewed as public edges, allowable to all players, and supplemented with a certain number of zero-cost edges allowable to only one player. These private edges are added and connected to the public ones in the following manner. First, a single-route network is created for every strategy of every player by connecting all the public edges that uses only in that strategy (that is, resources of the second and third kinds considered above) alternately with private edges for that player, such that the first and last edges in the route are private. (The order of the public edges is immaterial.) Second, the resulting single-route networks are connected in parallel, and the resulting network is connected in series with a similar single-route network in which all the public edges that uses in all of his strategies (that is, resources of the first kind above) alternate with private edges. After this is done for all players, the origin vertices of the networks are identified and so are their destination vertices. Clearly, for each strategy of each player in the congestion game defined above, all the resources that are included in that strategy, and only them, are also included in one of 's allowable routes in the network. Hence, the unweighted network congestion game is trivially isomorphic to the congestion game, and hence to .

Another variant of network congestion games capable of representing all finite games is weighted network congestion games (in the wide sense) without self-effect [START_REF] Milchtaich | Network topology and equilibrium existence in weighted network congestion games[END_REF]. These games are similar to "normal" weighted network congestion games in that a single nondecreasing cost function is associated with each edge , which together with the flow determines the cost of for all users. However, the cost of for a user is . It thus depends only on the total weight of the other users of , and may consequentially be different for users of different weights. In this respect, weighted network congestion games without self-effect are similar to unweighted network congestion games with player-specific costs. Using arguments similar to those given in the proof of Theorem 1, it is not difficult to show that every finite game can be represented as such a game.