Yves Elskens
email: yves.elskens@univ-amu.fr

D F Escande

F Doveil

Vlasov equation and N-body dynamics How central is particle dynamics to our understanding of plasmas ?

Keywords: basic plasma physics, Debye shielding, Landau damping, Vlasov equation, Nbody dynamics, Liouville equation, statistical ensemble PACS numbers : 52.20.-j Elementary processes in plasmas 05.20.Dd Kinetic theory 45.50.-j Dynamics .35.Fp Plasma: electrostatic waves and oscillations

Difficulties in founding microscopically the Vlasov equation for Coulomb-interacting particles are recalled for both the statistical approach (BBGKY hierarchy and Liouville equation on phase space) and the dynamical approach (single empirical measure on one-particle (r, v)-space). The role of particle trajectories (characteristics) in the analysis of the partial differential Vlasov-Poisson system is stressed. Starting from many-body dynamics, a direct derivation of both Debye shielding and collective behaviour is sketched.

a direct validation of the Vlasov-Poisson system in the continuum limit -actually, it bypasses the Vlasov picture by leading also easily to the analysis of collisions, and shows moreover how Debye shielding alters the Poisson description of the interactions [START_REF] Escande | Uniform derivation of Coulomb collisional transport thanks to Debye shielding[END_REF].

In this paper, we discuss some of the respective merits of the partial differential equation viewpoint inherent to the Vlasov equation and of the ordinary differential equation viewpoint of many-body Coulomb force dynamics. We hope this discussion will help the reader to grasp how delicate foundational issues the Vlasov equation involves. While our selection of topics is bound to remain far from exhaustive, we hope it will refresh the reader's view of the Vlasov equation. Specifically, Sect. II focuses on the meaning of the entity obeying the Vlasov equation, namely the distribution function. Section III comments on the role of particle trajectories in the theory of Vlasov-Poisson systems. Section IV sketches a specific instance where attention to trajectories provides better insight into well-known vlasovian phenomena -this is the only section where a specific microscopic model (actually, jellium) is considered.

II. STATUS OF THE DISTRIBUTION FUNCTION AND ITS EVOLUTION

Major merits of the Vlasov equation include 1. a rather successful application to the modeling of hot plasmas, 2. enabling physicists to use powerful [45] techniques of partial differential equations to address physical phenomena, 3. contributions to progress in functional analysis [46] motivated by challenges it raises, 4. stressing an interesting stage (in order to formulate a physical problem, one needs not merely to write its equations but first to choose the mathematical entities called to play) for describing the physics.

Issue 1 indicates how much needed is a comprehensive foundation for this fruitful model. We address issue 2 (and indirectly issue 3) in Sect. III. Let us stress in the current section the conceptually crucial issue 4 : the central quantity in the classical theory of the Vlasov equation is a non-negative continuous function f on Boltzmann's µ-space [START_REF] Penrose | Foundations of statistical mechanics -A deductive treatment[END_REF], f (r, v, t), which in the simple Coulomb interaction case obeys a closed evolution equation [47].

A. Pedestrian approach Classical, continuously differentiable distribution functions f (r, v, t) are often interpreted in terms of some averages. Since the particles are completely described by their positions and velocities, (r j (t), v j (t)), one first associates with these data a "spiky" distribution, involving Dirac distributions on (r, v) space. This is too wild an object for simple calculations, and one would gladly smoothe it. Moreover, as the field acting on particles varies in space, distinct particles are subject to different accelerations, and their motions may significantly separate with time. Yet, while particles drift apart from each other, other particles may come around them in such a way that the overall distribution does not seem to change much : recall a cloud made of many tiny drops, or the large-scale hydrodynamics of air made of many molecules.

Simultaneously, the field generated by the particles varies wildly on microscopic scales, so that one expects particle motions to be hardly accessible to the theory on these scales. Therefore, one would also gladly limit the particle motion analysis to large scales, over which the Coulomb field might appear smoother.

A first smoothing procedure is simply to replace the Coulomb interaction (with point sources) by a regularized, or mollified [START_REF] Neunzert | An introduction to the nonlinear Boltzmann-Vlasov equation, Kinetic theories and the Boltzmann equation[END_REF], interaction, where the source of the Coulomb field is a sphere with finite charge density centered on the point particle. The force on a particle is then also computed by summing the electric field over the sphere. This is done in particle-based numerical schemes, and if one keeps a fixed regularizing form function, one can even obtain the Vlasov equation from the N -body model rigorously (see e.g. [START_REF] Neunzert | An introduction to the nonlinear Boltzmann-Vlasov equation, Kinetic theories and the Boltzmann equation[END_REF]). However, the derivation of the Vlasov equation depends crucially on the size of the mollifying sphere. [48] A second procedure to formalize the "cloud" picture is to deem irrelevant some subtleties of particle motion. Rather, one pays attention to the evolution of the particles currently in a "mesoscopic" domain ∆U r (with size |∆U r | = ∆x∆y∆z), with velocities in a similar range ∆U v (with

|∆U v | = ∆v x ∆v y ∆v z ). The distribution function is then used to compute the "coarse-grained" distribution |∆U r × ∆U v | -1 ∆Ur×∆Uv f (r, v, t) d 3 r d 3 v
, where the range and domain of interest are large enough to contain so many particles (say Q ≫ 1) that their number would fluctuate moderately with time (say on the smaller scale Q 1/2 if these fluctuations follow a central-limit type of scaling, or Q 2/3 for a surface-vs-volume scaling in position space, or Q 5/6 for a boundary-vs-volume scaling in (r, v)-space).

To extract a smooth function f from this fluctuating picture, a third procedure consists in introducing an "ensemble" of realizations of the plasma, see e.g. [START_REF] Akhiezer | Non-linear theory and fluctuations[END_REF][START_REF] Hazeltine | The framework of plasma physics[END_REF][START_REF] Ichimaru | Basic principles of plasma physics -A statistical approach[END_REF]. The function f in the Vlasov equation is then viewed as the average of individual spiky distributions. As the force acting on a particle is due to the field generated by all other particles, this force is expressed as an integral over the distribution of those source particles, from which the target particle is excluded. This leads to computing a field E(r 1 , t) from an integral over the two-particle joint distribution f (2) (r 1 , v 1 , r 2 , v 2 , t), in the spirit of derivations of the Boltzmann equation in gas theory.

[49] Because of the long range nature of the Coulomb interaction, one expects particles to be "almost" independent and f (2) to almost factorize -otherwise, one should solve an evolution equation for f (2) where the source involves f (3) , etc... The resulting set of equations, viz. the BBGKY hierarchy, stands at the core of the statistical approach below. [50] Let us stress again that the usual conceptual setting for these derivations involves probabilistic averages over ensembles to generate smooth functions. Yet a physical plasma is a single realisation of the possible plasmas considered in an ensemble. The particles in it do not respond to the average field generated by the ensemble, and each particle follows a single, regular enough, trajectory (which is not a diffusion process).

An additional difficulty met with continuous distribution functions f solving the Vlasov equation is the absence of H-theorem. Indeed, the Vlasov equation preserves all functionals of the form R 6 G(f (r, v))d 3 rd 3 v, and evolutions towards a kind of equilibrium can only lead to the formation of finer and finer filaments rippling the surface representing f over (r, v) space.[51] When such ripples become finer than a typical interparticle distance in the N -body system, they lose physical significance [START_REF] Firpo | Long-time discrete particle effects versus kinetic theory in the self-consistent single-wave model[END_REF]. However, the BBGKY evolution equations do not incorporate such a destruction of unphysical ripples.

In numerical simulations, this filamentation is a delicate issue. On the one hand, modeling accurately the partial differential equation requires increasing computational power as filamentation proceeds. On the other hand, the numerical smoothing due to various interpolation schemes is not granted to reproduce the physical smoothing of the distribution function due to phenomena not included in the smooth Vlasov model, such as finite-N effects, perturbing interactions, etc.

To conclude, there is a single Vlasov equation, but there are various views of the distribution function whose evolution it is meant to describe. For a given plasma physics problem, which of these views, if any, should be considered ? This issue is usually overlooked, and the outcome of the vlasovian calculation is deemed relevant.

B. Technical approach

The statistical approach to deriving the Vlasov equation (see e.g. Appendix A in [START_REF] Ichimaru | Basic principles of plasma physics -A statistical approach[END_REF]) starts from N -body dynamics, introduces the high-dimensional phase space Γ = R 6N , and considers the Liouville equation for a statistical measure f on Γ. [52] This statistical measure may be interpreted as the probability distribution of, say, N "replicas" (N ≫ 1) of the N -body system, and the symmetrized one-particle marginal[53] of f obeys the first equation in the BBGKY hierarchy.

One then identifies this symmetrized marginal with a one-particle measure on the "molecular" µ-space[54], and one expresses the force field generating the evolution by some integrals of the two-particle marginal of f. For the Vlasov-Poisson system, one then requires that the electric field in the Vlasov equation solves the Poisson equation where the source is the average of the individual particle distributions generated by all replicas.[55] In this approach, the evolution of f is thus subsidiary to the evolution of many replicas, each of which evolves independently of the N -1 other replicas, and the f of interest in the Vlasov equation is then an average over N replicas in the limit N → ∞. But how can (N -1) thought-experimental replicas drive the evolution of the single physically realized N -body system ? How should the field acting on a given particle in the physically observed plasma be the by-product of an average over N Gedanken plasmas ?

The dynamical approach to the Vlasov equation [START_REF] Kiessling | The microscopic foundation of Vlasov theory for jellium-like newtonian N -body systems[END_REF][START_REF] Spohn | Large scale dynamics of interacting particles[END_REF] starts with an actual system of N bodies, interacting by instant-action-at-a-distance or via dynamical fields (see e.g. [START_REF] Firpo | Kinetic limit of N -body description of wave-particle self-consistent interaction[END_REF] for a simple example). For undistinguishable particles, the data of N points in µ-space, namely a set[56]

M = {(r 1 , v 1 ), . . . , (r N , v N )}
(thanks to Coulomb repulsion, no particles can be at the same position, hence this set counts exactly N points), is equivalent to the counting measure dµ

(N ) = N j=1 δ(r -r j )δ(v -v j ) d 3 r d 3 v on R 6 .
When particles move with time, the measure also evolves, and for finite N the evolution of µ (N ) provides all information on the motion of all particles, as particles cannot swap their identity (exchanging two particle labels would require their trajectories to be discontinuous, or to meet at a same position with the same velocity -which is ruled out by the dynamics). The measure µ (N ) determines the force field exactly as the N particle data do, and this field generates the vector field according to which µ (N ) is transported. In order to handle the limit N → ∞, it is convenient to consider the normalized empirical measure N -1 µ (N ) , which is non-negative and verifies N -1 µ (N ) (R 6 ) = 1. Indeed, this makes µ (N ) /N formally akin to a probability measure, and indeed one may interpret N -1 µ (N ) (A) as the fraction of the plasma particles which are in some subset A ⊂ R 6 , or as the probability that a particle with randomly picked label[57] be in A.

The limit N → ∞ makes sense formally for non-negative normalized measures, and indeed this space of measures can be equipped with various kinds of distances generating physically reasonable topologies [START_REF] Mischler | A new approach to quantitative propagation of chaos for drift, diffusion and jump processes[END_REF][START_REF] Spohn | Large scale dynamics of interacting particles[END_REF] to give an operational meaning to the notation lim N →∞ . One such distance was considered by Kolmogorov and Smirnov to test the likelihood for sample random data to follow a given law, and it is used in convergence theorems for distributions [START_REF] Dudley | Real analysis and probability[END_REF]. On the contrary, the limit N → ∞ is ill-defined in the phase-space approach involved in the Liouville equation, on which the BBGKY hierarchy relies. The very space in which the phase point (representing the N -particle system) evolves varies with N : its dimension increases like N . Therefore, phase space Γ is simply not the good stage for performing the limit.

The derivation of the Vlasov equation in the (dynamical) measure approach is rather short (actually, shorter than the BBGKY-hierarchy based derivation) and both conceptually and physically clear, provided the interaction is not too singular (with the statistical approach providing no better derivation in the singular case either). Key ingredients are

• the existence and uniqueness of solutions to individual particles' equations of motion in a given, regular enough, force field, say E(r, t) (with e.g. Lipschitz regularity, i.e.

E(r, t) -E(r ′ , t) ≤ K(t) r -r ′ ),
• by duality [58] and because the single-particle dynamics in the given field is measurepreserving in the "molecular" µ-space, the existence and uniqueness of evolution of measures in the same given, regular enough, force field E(r, t),

• regularity of the force field generated by arbitrary measures (which fails for Coulomb field with point sources),

• a self-consistency argument to formulate the Vlasov equation coupled with the fields as a fixed point problem in a suitable (Banach) function space -much like the standard technique used for solving iteratively ordinary differential equations.

The measure approach reaches further than just proving that the limit N → ∞ commutes with time evolution. It actually shows that, given initial measures µ 1 (0) and µ 2 (0) on (r, v)space, which need not be empirical nor absolutely continuous, the Vlasov equation defines unique time-dependent measures µ 1 (t) and µ 2 (t) and that the distance between the evolved measures is controlled by their initial distance and a constant C depending on µ 1 (0) and µ 2 (0),

dist(µ 1 (t), µ 2 (t)) ≤ C ′ dist(µ 1 (0), µ 2 (0)) e C|t| . (1)
Here C ′ may occur for technical reasons, and C is an upper estimate for the largest Lyapunov exponent, controlling the rate at which particles separate in the fields E 1 and E 2 generated by the matter distributions described by µ 1 and µ 2 (see [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF][START_REF] Firpo | Kinetic limit of N -body description of wave-particle self-consistent interaction[END_REF] for a smooth interaction example). This estimate implies that, to ensure dist(µ (N ) (t)/N, µ f (t)) ∼ δ with dist(µ (N ) (0)/N, µ f (0)) ∼ N -c (with c = 1/2 or 1, say), one needs an initial accuracy corresponding to some N ∼ (C ′ /δ) 1/c exp(C|t|/c), which is too demanding when the times of interest exceed a few Lyapunov e-folding times.

As physicists are more interested in estimating specific observables, say local current densities or electric field energy density, estimates like (1) provide upper bounds on the errors due to approximating µ 1 (0) with µ 2 (0). However, actual discordances between the values of specific observables are generally much smaller -for instance, if µ 1 (0) and µ 2 (0) describe two stationary solutions, the observables are constants of the motion.

III. VLASOV-POISSON EVOLUTION

While the dynamical foundation of the Vlasov equation, in the previous section, stressed the motion of particles, we must stress that some important results on the Vlasov-Poisson system make no use of particle trajectories. In this section, we first indicate how the smooth solutions to the Vlasov-Poisson system may involve non-trajectorial concepts, and then comment on some aspects of the system which may, or need to, involve trajectories.

A. Without trajectories

As a partial differential equations system, the Vlasov-Poisson model is well-posed and much studied (see e.g. [START_REF] Rendall | Kinetic equations[END_REF]). Theorems on solutions existence rely on some regularity in initial data, typically f has compact support and bounded, continuous derivatives ∂ r f , ∂ v f . Proofs involve estimates for the spatial density ρ(r, t) = R 3 f (r, v, t) d 3 v in various L p norms, along with estimates for the electric field E (whose L 2 norm is proportional to the total potential energy) and its gradients, which involve norms of ∂ r ρ. Moments of f (and of course of ρ), as well as the upper bound on particle velocity P (t) = sup{|v| : f (r, v, t) > 0 for some r}, play a role in the analysis of the Vlasov-Poisson system.

It may occur that gradients of f get steep, due to filamentation in µ-space. Yet, gradients of ρ or of f (v, t) := R 3 f (r, v, t) d 3 r may be better controlled, and even decay in some sense to zero, as e.g. in the nonlinear theory of Landau damping [START_REF] Bedrossian | Landau damping : paraproducts and Gevrey regularity[END_REF][START_REF] Mouhot | On Landau damping[END_REF]35,[START_REF] Villani | Particle systems and nonlinear Landau damping[END_REF].

The analysis of continuous solutions to the Vlasov-Poisson system does not rely only on the more usual physical quantities like ρ or j(r, t) = R 3 vf (r, t) d 3 v, and on standard global invariants like total momentum or energy. It also takes advantage of Casimir invariants, which may be expressed in terms of the functions

C a [f ](t) := R 6 1(f (r, v, t) > a) d 3 r d 3 v (2)
where 1(A) = 1 if A is true and 1(A) = 0 otherwise. Note that C a is a decreasing function of a, with C a = +∞ for any a < 0, and C 0 being the measure of the support of f . These functions C a measure the level sets of f and provide a tool for calculating other functionals like

f L 1 = ∞ 0 C a [f ] da, f p L p = p ∞ 0 a p-1 C a [f ] da, or R 6 f (r, v) ln(f (r, v)/c)d 3 rd 3 v = ∞ 0 (1 + ln(a/c))C a [f ]da.
The Vlasov equation preserves all Casimir invariants, and a physical interpretation for these conservation laws is that vlasovian evolution must be invariant under the group of particle relabelings. Such a group is discrete for the N -body system (and then C a is undefined), but it is continous for smooth f , and one expects it to generate integral invariants following Noether's theorem. Adding to the total energy a suitably chosen Casimir invariant enabled proving the stability of some spherically symmetric equilibria of the gravitational Vlasov-Poisson system [START_REF] Mouhot | Stabilité orbitale pour le système de Vlasov-Poisson gravitationnel[END_REF].

In numerical simulations and in simple analytic models, the case where C a is a simple step function occurs for waterbag distributions [START_REF] Bertrand | Non linear electron plasma oscillation : the "water bag model[END_REF], equal to h on a domain A(t) with measure 1/h, and vanishing outside A(t) (then C a = 0 for a ≥ h and C a = 1/h for 0 ≤ a < h).

B. With trajectories

Particle trajectories are instrumental in understanding the Vlasov equation because this equation transports the distribution f along the characteristics of the Vlasov operator. Good control on particle trajectories is crucial e.g. to proofs of the existence of solutions globally in time [START_REF] Horst | On the classical solutions of the initial value problem for the unmodified non-linear Vlasov equation -I. General theory[END_REF][START_REF] Pfaffelmoser | Global classical solutions of the Vlasov-Poisson system in three dimensions with generic initial data[END_REF][START_REF] Schaeffer | Global existence of smooth solutions to the Vlasov Poisson system in three dimensions[END_REF] (using the Lipschitz regularity of the electric field generated by the distribution function) and to the proof of Landau damping in the nonlinear regime [START_REF] Bedrossian | Landau damping : paraproducts and Gevrey regularity[END_REF][START_REF] Mouhot | On Landau damping[END_REF]35,[START_REF] Villani | Particle systems and nonlinear Landau damping[END_REF].

Estimates on individual particle trajectories in small field limits are obtained for the fields generated by the Vlasov-Poisson dynamics, and provide further insight in the latter, using estimates for velocity moments of the distribution function [START_REF] Lions | Propagation of moments and regularity for the 3-dimensional Vlasov-Poisson system[END_REF][START_REF] Ch | About the velocities of self-interacting particles in the 1-D torus[END_REF]. Estimates on particle trajectory crossings are also used to construct Lyapunov functionals for the evolution from "small" initial data in order to assess the decay to spatial uniformity and zero electric field in long times, and the "asymptotic completeness"[59] of dynamics [START_REF] Ha | New Lyapunov functionals of the Vlasov-Poisson system[END_REF][START_REF] Choi | Asymptotic behavior of the nonlinear Vlasov equation with a self-consistent force[END_REF].

Incidentally, Mouhot and Villani also stress the stabilizing role of plasma echoes (a particle effect, associated with bunching) for the nonlinear theory of Landau damping. However, characteristics (i.e. trajectories) are not considered individually ; rather, smoothness of f is important, and nonlinear particle motion does not break it for finite time -while mixing in (r, v) space generates small scale oscillations (ripples) which are homogenized in the long run.

On the other hand, particle trajectories may be ignored when considering stationary solutions, see e.g. the orbital stability problem solved by Lemou, Méhats and Raphaël (see [START_REF] Mouhot | Stabilité orbitale pour le système de Vlasov-Poisson gravitationnel[END_REF] for a review).

Nor do they appear explicitly in the discussion of Bernstein-Greene-Kruskal modes and trapping in space-periodic Vlasov-Poisson system [START_REF] Zh | Small BGK waves and nonlinear Landau damping[END_REF].

Note that all these statements refer to the Vlasov-Poisson system as a partial-differential deterministic evolution equation system, in which any particle is following a trajectory defined by the Vlasov equation characteristics. It is not analysed with probabilistic methods.

IV. DIRECT MANY-BODY DYNAMICS

A finite-N approach always stresses trajectories. Analogues of the stability theorems for Landau damping in the nonlinear regime would be a Kolmogorov-Arnol'd-Moser theorem [START_REF] Bedrossian | Landau damping : paraproducts and Gevrey regularity[END_REF] for globalin-time stability, or a Nekhoroshev-type theorem for stability over exponentially long times (say exp(aN b ) for some a, b > 0). This is indeed a challenging target for mathematical analysis (see the recent tutorial [START_REF] Villani | Particle systems and nonlinear Landau damping[END_REF]).

Interestingly, (partial-differential-equations-based) Landau damping theorems involve periodic boundary conditions, so that the electrons may be viewed as immersed in a neutralizing background (the "jellium" model) -while results on asymptotic completeness involve small initial data in an infinite space.

For periodic jellium, with finitely many bodies, interacting through repulsive Coulomb force, no two particles can be simultaneously at the same point. The equations of motion are well-posed, and solutions to the initial value problem (Cauchy problem) exist for all times.[60] The continuum limit, where N → ∞ along with µ (N ) /N approaching an absolutely continuous measure with a smooth density f , is however difficult because the limit imposes that some particles come arbitrarily close to each other, and the characteristic time scales for their relative motion become arbitrarily short.

The short-range singularity of the Coulomb potential is too strong for the usual rigorous derivations of the Vlasov equation to apply. Kiessling [START_REF] Kiessling | The microscopic foundation of Vlasov theory for jellium-like newtonian N -body systems[END_REF] exposes a new approach along with a review of the state of the art.

Without relying on the Vlasov description, one can use the N -body description directly to understand some fundamental physical effect in plasmas. This issue was recently revisited [START_REF] Escande | New foundations and unification of basic plasma physics by means of classical mechanics[END_REF] for N electrons (with mass m e and charge -e) in a cube with side L with periodic boundary conditions (the latter play the role of a neutralizing background). For large N , when the particle positions are not too far from a spatially uniform distribution, and velocities are close to a distribution f 0 (v), one separates the motion of a particle into a ballistic approximation, r (0) l (t) = r l0 + v l0 t, and a correction δr j (t) = r j (t)r (0) l (t) due to the interactions. Perturbative analysis [61] shows that the dominant contribution to their interaction is not given by the full Coulomb potential (with Fourier

components φ(m, t) = -N j=1 e/(ǫ 0 k m 2 ) exp(-ik m •r j (t)) for k m = 2πm/L with m ∈ Z 3 \{0}),
but only the Debye-shielded potential δΦ(r, t) = j δΦ(rr j (0)ṙj (0)t, ṙj (0)), where

δΦ(r, v) = - e L 3 ǫ 0 m =0 exp(ik m • r) k m 2 ǫ(m, k m • v + iε) (3)
and ǫ approaches[62] the usual dielectric function

ǫ(m, ω) = 1 -e 2 /(L 3 m e ǫ 0 ) (ω -k m • v) -2 f 0 (v) d 3 v (4)
in the N → ∞ limit.

While Debye shielding of probes or test charges is well known in the kinetic theory of plasmas, this result brings forward the direct appearance of Debye shielded Coulomb potential in the dynamics -which is not usually stressed in classical Vlasov-Poisson theory. Moreover, the faster decay of the shielded force enables one to compute transport phenomena in terms of particle "binary collisions" over the whole range of impact parameters [START_REF] Escande | Uniform derivation of Coulomb collisional transport thanks to Debye shielding[END_REF], while the classical collisional models were constructed separately for large and small impact parameters. In the same calculations, the collective response of particles to any one of them also "dresses" the latter, so that the effective description of the plasma deviates simultaneously from both the Vlasov equation (stressing now N bodies) and from the Poisson equation (through shielding and dressing).

The potential (3) lends itself easily to Laplace analysis in time [START_REF] Escande | New foundations and unification of basic plasma physics by means of classical mechanics[END_REF]. In the N → ∞ limit, it then also approaches formally (outside its sources) the potential given by the standard expression involving the Landau contour calculation, for both the unstable (growing wave) case and the damping case. The twist in this derivation of Landau damping in the linear regime is that it involves not the full Coulomb interaction but its (perturbatively constructed) shielded version.

This comes as a surprise, for screening generates a "short range" effective model ; yet the shielded potential is still Coulomb-singular at its source, so that the standard derivations of the Vlasov equation do not apply better than for the full Coulomb interaction.

V. FINAL REMARKS

Let us first recall that finite-N dynamics may always depart from smooth Vlasov solutions on the characteristic time scales of instabilities, as noted in eq. ( 1). Apart form this obvious example, it is known that for interactions smoother than Coulomb's, finite-N distributions can indeed remain quasi-stationary for extremely long times (see [START_REF] Levin | Nonequilibrium statistical mechanics of systems with long-range interactions[END_REF] for a recent review), even if the finite-N model does not admit corresponding periodic solutions [START_REF] Elskens | Finite-N dynamics admit no traveling-wave solutions for the hamiltonian XY model and single-wave collisionless plasma model[END_REF]. Moreover, the large, finite N regime of near-equilibrium distributions, as occurs in the physical problem motivating research on orbital stability, is of interest for the Newton gravitation interaction [START_REF] Mouhot | Stabilité orbitale pour le système de Vlasov-Poisson gravitationnel[END_REF], where it also branches to the elegant special solutions to the N -body dynamics with symmetries known as "choreographies" [START_REF] Montgomery | N-body choreographies[END_REF]. Similar results in the Coulomb case relate to the motion of cold ion clouds or non-neutral plasmas [START_REF] Dubin | Trapped nonneutral plasmas, liquids, and crystals (the thermal equilibrium states)[END_REF][START_REF] Kiessling | A note on classical groound state energies[END_REF].

As physicists have since long understood that finite-N systems are never completely described by classical Vlasov-Poisson theory, long-time evolutions of plasmas are often described using other kinetic models in the N → ∞ limit. Rigorous derivations of such models, like the Landau equation or the Balescu-Lenard-Guernsey equation, are a challenge for this century (see [START_REF] Lancellotti | From Vlasov fluctuations to the BGL kinetic equation[END_REF] and the truly short ch. 6 in [START_REF] Spohn | Large scale dynamics of interacting particles[END_REF]). Actually, corrections to the Vlasov model often go by the name of "collisions", and we noted in Sect. IV how surprising they may be even in their simpler instances [START_REF] Escande | New foundations and unification of basic plasma physics by means of classical mechanics[END_REF][START_REF] Escande | Uniform derivation of Coulomb collisional transport thanks to Debye shielding[END_REF].

Collisions do not necessarily account for the graininess effects due to finite N [START_REF] Firpo | Long-time discrete particle effects versus kinetic theory in the self-consistent single-wave model[END_REF]. [46] This progress makes the modern theory of partial differential equations incommensurately more powerful, benefiting to many other fields.

[47] The Poisson equation is not an evolution equation : it simply provides an explicit relation expressing the electric field in terms of the current, instantaneous f -in this sense, the electric field does not describe additional degrees of freedom, but is a mere subsidiary entity slaved to the particles.

[48] Physically, one expects the classical model to fail anyway at describing Coulomb-interacting particles at very short distances. Indeed, the "classical radius of the electron" (at which scale the point particle model should break down) is far smaller than its Compton length, so that quantum effects are expected to imply a different modeling. It may thus be considered physically sensible to model the N -body plasma with a "quantum-mollified" Coulomb interaction (even for large, finite N ). This will seem even more harmless since the quantum description of two-body interaction [START_REF] Bohm | Quantum theory[END_REF] leads to a scattering cross section identical to the classical Rutherford cross section if particles bear charges with equal signs -moreover, the same expression applies to particles with opposite signs, for which the classical picture allows for arbitrarily close approaches. Besides, physicists most often consider plasmas where typical binary collisions involve closest approaches much further away than the Compton length, and molecular dynamics simulations simply integrate the Coulomb force for short range while smoothing the interaction for longer range [9]. Yet the mathematical issue is important, for it should provide insight into more complicated problems too.

[49] Until 2013, rigorous derivations of the Boltzmann equation were limited to short times (on the order of the mean free path) or small initial data (ultimately expanding in vacuum without collision). This is a topic of hard and active research, see [START_REF] Gallagher | From Newton to Boltzmann : hard spheres and short-range potentials[END_REF][START_REF] Pulvirenti | On the validity of the Boltzmann equation for short range potentials[END_REF] and references therein.

[50] Note that, in the N -body case, giving the spiky distribution amounts to locating all particle, and the joint distribution for pairs is uniquely determined by these data. In contrast, given a smooth distribution f (1) on R 6 , there exist infinitely many joint distributions f (2) on R 12 such that f (2) (r 1 , v 1 , r 2 , v 2 , t) d 3 r 2 d 3 v 2 = f (1) (r 1 , v 1 , t) along with f (2) (1, 2) = f (2) (2, 1). The factorization f (2) (1, 2) = f (1) (1)f (1) (2) may hold at initial time but is generally not preserved by the time evolution. This "propagation of (initial) chaos" (or "Kac property") is an important, non-trivial issue, for understanding the N → ∞ limit [START_REF] Mischler | A new approach to quantitative propagation of chaos for drift, diffusion and jump processes[END_REF].

[51] This is also viewed from the fact that the Vlasov equation for hamiltonian dynamics "shuffles" (without altering their volumes) the respective domains over which f takes its various values, i.e. at any time t the set {(r, v) : a < f (r, v, t) ≤ b} has the same 6-dimensional volume (see [START_REF] Hazeltine | The framework of plasma physics[END_REF] for this shuffling image).

For the one-dimensional models, with (r, v) = (x, v), volume preservation reduces to area preservation, which makes hamiltonian dynamics so different from dissipative ones. Recall also how mixing occurs in non-diffusive, viscosity-free incompressible fluid flows.

[52] Given an initial data z ∈ Γ, the N -body equation of motion ż = g(z) generates the evolution of the 6N degrees of freedom so that z(t) = T t,0 (z 0 ) for initial data z(0) = z 0 . A statistical measure

It is a pleasure

  for YE to acknowledge stimulating discussions with G. Belmont, F. Califano, C. Krafft, G. Manfredi, Ph. Morrison, F. Pegoraro, F. Valentini and participants to Vlasovia 2013 in Nancy. The authors are grateful to L. Couëdel for a critical reading of the manuscript, and to their colleagues in Marseilles for fruitful comments. of variables.

on Γ is a distribution f, and the Liouville equation ∂ t f + ∂ z • (gf) = 0 transports an initial measure f 0 for z 0 into an evolved measure f t = f 0 • T 0,t (just by finding the preimage at time 0 of the state at time t). For instance, for free motion, z = (r

integration with the Dirac distribution generates the one-particle distribution for particle j, and the sum symmetrizes over all particles.

[54] This space would be the phase space of a particle, if the force field acting on it were known (though possibly time-dependent). However, for the many-body problem, particles interact, so that the phase space is Γ = ⊗ N j=1 R 6 while µ-space is just R 6 (see notes 72 and 118 in [START_REF] Ehrenfest | The conceptual foundations of the statistical approach in mechanics[END_REF] for a historical gas theory analogue).

[55] The source of E is not the N -1 N j=1 δ(rr j ) given by a single plasma realisation, but the density

obtained from an average over N replicas.

[56] Not a sequence where labels do matter !

[57] One shall not confuse this "label" randomness, essentially related to particle undistinguishability in a single realization of the plasma, with the random choice of a replica out of N possible realizations of the plasma.

[58] This is just using the fact that µ(T t,0 (y) ∈ A) = µ(y ∈ T 0,t A), because the space of measures µ is dual to the space of observables q, in the linear algebraic interpretation of the formula R 6 q(y) dµ(y) [START_REF] Appel W (transl | Mathematics for physics and physicists[END_REF].

[59] This issue is whether the nonlinear Vlasov flow can or cannot be approximated by the free streaming flow in the limit t → ∞.

[60] For actual plasmas, with ions and electrons, Coulomb attraction permits that particles come arbitrarily close, with potential energy diverging to -∞ while kinetic energy diverges to +∞. This makes their N -body description even more difficult. The similar difficulty in celestial mechanics is well known [START_REF] Falcolini | Collisions and singularities in the n-body problem[END_REF].

[61] See appendix B in [START_REF] Escande | New foundations and unification of basic plasma physics by means of classical mechanics[END_REF] for the outline of a nonlinear version.

[62] The iε prescription stems from inverting the Laplace transform as the integral in ǫ is undefined for real-valued ω = k m • v. As this diverging ǫ occurs in the denominator in (3), however, the finite N version does not need the iε. A similar discussion occurs in the analysis of discrete analogues of van Kampen-Case modes [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF].


Benjamin Dollet

Isabelle Cantat

Deformation of soap films pushed through tubes at high velocity

 holds up to surprisingly high velocities, at which the capillary and Weber numbers are no longer small parameters. The tube geometry is varied and accounted for through the notion of hydraulic diametre. A limit of stability of the film, beyond which the films burst or evolve unsteadily, is predicted, and captures quantitatively the observations. The new questions raised by our results on the dissipation in soap films are discussed, especially the role of Plateau borders and inertial effects.

Introduction

Aqueous foams are widely used in industry, in areas as diverse as ore flotation, oil extraction, food and pharmaceutical production and blast noise reduction, to cite a few.

Many foam flows of practical interest are dominated by viscous friction because of high velocity gradients, either in bulk, or localised close to confining boundaries. For instance, foams in porous media, used for enhanced oil recovery [START_REF] Hirasaki | Mechanisms of foam media: Apparent viscosity in smooth capillaries[END_REF] and soil remediation [START_REF] Chowdiah | Foam propagation through soils for enhanced in-situ remediation[END_REF], often flow through narrow gaps as a train of soap films, in contact with the confining walls through Plateau borders (PBs) where viscous dissipation takes place and dictates the foam flow properties [START_REF] Rossen | Percolation theory of creation and mobilization of foams in porous media[END_REF][START_REF] Kornev | Foam in porous media: thermodynamic and hydrodynamic peculiarities[END_REF]. The influence of confinement is also paramount in microfluidics, where foams are used in lab-on-a-chip applications (see [START_REF] Marmottant | Microfluidics with foams[END_REF] for a review).

Foam flows dominated by viscous friction are also of fundamental interest: foam rheology has received a considerable attention in recent years [START_REF] Höhler | Rheology of liquid foams[END_REF], but mostly at velocities small enough for the foam to remain in equilibrium configurations, i.e. for structural equilibrium rules to apply. However, deviations from these equilibrium rules have been recently evidenced [START_REF] Drenckhan | Rheology of ordered foams-on the way to discrete microfluidics[END_REF][START_REF] Besson | Dissipation in a sheared foam: From bubble adhesion to foam rheology[END_REF], due to bubble/bubble or bubble/wall viscous friction. This calls for more studies of foam rheology at high velocity gradients, where the structure is brought out of equilibrium because of a high viscous resistance.

In this paper, we study the most elementary experiment of high-velocity foam rheology in confined media: the fast flow of an isolated soap film pushed through a tube of constant section. We especially quantify over a full range its shape deviation from the simple equilibrium configuration, perpendicular to the confining walls. We extract a precise experimental criterion for its rupture in terms of critical velocity.

Materials and methods

A glass tube was connected through flexible plastic tubes to a pressurised nitrogen bottle, which supplies nitrogen through an electronic flow-rate controller (Brooks) (Fig. 1).

We have used three circular glass tubes, of inner diameters 3.6, 5.9 and 8.8 mm, and a square one, of side 10.0 mm with rounded edges of radius of curvature lower than 1 mm.

The tube is tilted upwards with an angle of about 20 • with respect to the horizontal, so that the plastic tube to which it is connected is at lower level and can be used as soap reservoir. Soap solution is injected by a syringe, directly through the open end of the glass tube. When turned on, the gas flows into this reservoir, and creates soap films which are subsequently pushed through the glass tube. Transient effects were also investigated using films prepared first at low velocity, and then accelerated (Sec. 4.2). When starting an experiment, we waited for the solution to wet entirely the glass walls, because of the first soap injections and films. The wetting film remains intact during the whole experiment: the signature of dewetting would be the nucleation and growth of dry patches and an irregular behaviour of the films that was never observed.

The soap solution is a solution of SDS (Sigma-Aldrich) dissolved in a mixture of ultrapure water (Millipore) and of 10% wt glycerol (Sigma-Aldrich). The concentration of SDS is 10 g/l, above the critical micellar concentration (cmc) of 2.3 g/l. Its bulk viscosity, measured in a Couette rheometer (Anton Paar) is η = 1.2 mPa s. Its surface static and dynamic properties were measured with a tensiometer (Teclis) by the rising bubble and oscillating bubble method; the surface tension is σ = 36.8 ± 0.3 mN/m and the surface modulus E S , defined as in [START_REF] Denkov | Wall slip and viscous dissipation in sheared foams: Effect of surface mobility[END_REF], was below the noise level of 1 mN/m at frequency 0.2 Hz. All solutions were used within a day from fabrication.

The motion of a soap film is recorded with a high-speed camera (Photron APX-RS) in the central zone of the glass tube, to avoid end effects. The main results, presented in Sec. 4.1, have been obtained with a 2 cm long field of view (measured along the tube axis). We checked that no significant variations of shape and velocity occur along this trajectory: hence, the data in Sec. 4.1 concern steady film motions. Additionnal experiments, discussed in Sec. 4.2, have been made on unsteady film motions, with a 15 cm long field of view to investigate the long range stability of the films, especially for the highest velocity values. Exemples of imaged soap films are displayed in Figs. 2 (circular tube) and 3 (square tube). From each movie, we plot as a spatiotemporal diagram (z, t), with z the streamwise axis and t the time, the evolution of the gray levels along the tube centreline (Fig. 4a). The soap film appears as a dark band, which slope α yields the film velocity v: v = dz/dt = cot α. The slope can be extracted with an accuracy better than δα = 0.1 • , which gives the uncertainty on v: δv = δα/ sin 2 α.

Above a given velocity of order 40 cm/s, the soap film shows a significant curvature pointing downstream (Figs. 2b,c and4b). To measure this curvature, we have performed the following home-made contour detection procedure. First, we have thresholded and binarised the images (Fig. 4c), then we have extracted the coordinates (ξ, z) of the points located on the curved profile (ξ = 0 at the centre of the tube). This set of coordinates was then fitted to the polynomial approximation of degree four of a circle of radius R:

-z = (ξ -ξ 0 ) 2 /2R + (ξ -ξ 0 ) 4 /8R
3 and we extract R as the best fit parameter (Fig. 4d).

The offset ξ 0 is required because, in some experiments, the soap film is not symmetric, as in Fig. 2b where the top part of the film lags behind. All this procedure assumes a constant curvature; we confirmed that this was a realistic assumption by checking the close agreement between the whole experimental trace and the best fit, which also gives a small relative uncertainty on the fitting parameter R (this uncertainty is represented by the vertical error bars in Fig. 5). Notice that the film curvature is not measurable at lower velocities, because the curved profile is hidden by the thicker trace of the PB (Fig. 2a).

Theory

Film deformation

We consider the steady motion of a soap film pushed at imposed velocity v through a cylindrical tube of arbitrary shape, with a perimeter P and a cross-section S. The film displays a curvature H because of the pressure drop ∆P associated to its motion, due to the viscous friction in the liquid. After Laplace law, the two latter quantities are related through: ∆P = 2σH, accounting for the fact that the soap films is bounded by two gas/liquid interfaces. Hence, if the pressure on both sides of the film are equilibrated, the film has a constant curvature. Here, we assume that the surface tension is constant in the transverse soap film (it needs not be in the PB; see Sec. 5.2). Following the assumption of [START_REF] Cantat | Dissipation in foam flowing through narrow channels[END_REF], we assume that the friction force of an elementary length dℓ of PB, of normal unit vector pointing downstream n ′ (Fig. 3), is proportional to its projected length across the tube, dℓ n ′ • e z , and to the friction force f (v) per unit PB length; therefore, the total friction force experienced by the soap film is Pf (v). Hence, the steady motion of the soap film obeys the equation: 2σSH = Pf (v).

In our experiments, close to the centre of the tubes, the film shows a constant radius of curvature R; therefore, H = 2/R. Let further b = 4S/P be the hydraulic diametre of the tube; for a circular tube of radius a, b = 2a, and for a square tube of side c, b = c. Thus, we have the simple prediction that R/b should be independent from the tube geometry:

R b = σ f (v)
.

(3.1)

Limit of steady film deformation

The above discussion assumed the existence of a steady shape of the soap film into motion. This clearly fails at high velocity, because the curvature of the film becomes too high for the film to meet the tube.

This idea is easy to quantify for a circular tube of radius a. A transverse film in steady motion meets then the wall tube with a constant angle ϕ, obeying sin ϕ = a/R (Fig. 4b), and f (v) = 2σ sin ϕ after (3.1). Hence, the minimal radius of curvature is simply: with a wall, say at y = 0, of length dℓ and of projected length dx (Fig. 3). This portion of PB meets the upstream and downstream wetting films, which tensions cancel out, and the transverse film, which two gas/liquid interfaces exert a tension in its locally tangential direction t along the length dℓ. Furthermore, it experiences a viscous force -f (v) dx e z .

R circle min = a = b 2 , ( 3
Hence, the force balance along e z on the portion of PB writes:

-f v dx + 2σ dℓ t • e z = 0. Moreover, t = n ∧ t ′ ,
2 x )(1 + ζ 2 x + ζ 2 y ). Given that dx/dℓ = 1/ 1 + ζ 2
x , the force balance gives the condition:

f (v) 2σ = ζ y 1 + ζ 2 x + ζ 2 y = cos γ, (3.3)
which is analogous to a capillary surface meeting a wall at prescribed contact angle γ = π/2 -ϕ, studied by [START_REF] Concus | On the behavior of a capillary surface in a wedge[END_REF]. They predicted that in a wedge of interior angle 2δ, there exists no bounded surface if γ + δ < π/2. In our case, the square tube has four corners with an angle δ = π/4, hence the condition of existence of a shape of a transverse film in steady motion writes simply: ϕ ≤ π/4 or, after (3.3), f (v)/σ ≤ √ 2.

Hence, after (3.1), the minimal radius of curvature is:

R square min = c √ 2 = b √ 2 . (3.4)
When the friction law f (v) is known, inserting R min in (3.1) yields a prediction of the maximal velocity of a stable soap film.

Friction law

Several investigations have studied the viscous friction associated to the motion of a soap film. In a seminal paper, [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] has proposed a theoretical expression for the friction associated to the motion of an infinitely long bubble in a tube, in the limit where inertia is negligible with respect to viscous and capillary effects, and where the capillary number Ca = ηv/σ, comparing viscous and capillary effects, is small. In that case, viscous friction comes from the transition zone matching the flat wetting film between the bubble and the wall, and the bubble tip. The length of the transition zone, and the thickness of the wetting film h ∞ , are slaved to the radius of the bubble tip, i. e.

the tube radius, but this cancels out in the expression of f (v):

f (v) = 4.70βσCa 2/3 , (3.5)
The factor β depends on the boundary condition at the gas/liquid interface: for a free shear boundary condition (mobile interfaces), β = 1, whereas for a no-slip boundary condition (rigid interfaces), β = 2 1/3 ≃ 1.26 (Appendix A). [START_REF] Hirasaki | Mechanisms of foam media: Apparent viscosity in smooth capillaries[END_REF] later extended this result to the case of soap films, showing that the relevant length becomes the PB radius, although this does not alter the expression of f (v). This idea was also used by [START_REF] Denkov | Wall slip and viscous dissipation in sheared foams: Effect of surface mobility[END_REF] We will compare our data to this formula in Sec. 4, before discussing its underlying limitations in Sec. 5, notably the possible role of the PB size and of inertial effects.

Experiments

Film deformation

We first plot for each circular tube R as a function on v, and we fit the data by a power law (inset of Fig. 5). We obtain exponents of -0.64 ± 0.08, -0.71 ± 0.04 and -0.71 ± 0.04 for tubes of respective radius 3.6, 5.9 and 8.8 mm, all compatible with -2/3. Next, fitting each series by a law A/v 2/3 , and plotting the best fit parameter A as a function of a for the three tubes, we have checked that the data follow a linear law A ∝ a within 6% (data not shown). Hence, we can plot R/2a as a function of v for all tubes (Fig. 5). The data collapse on a single curve, which best fit by a power law yields an exponent -0.65 ± 0.04.

Therefore, the scaling R ∝ a/v 2/3 predicted by (3.6) is experimentally verified, up to a velocity of 2.0 m/s, corresponding to a capillary number of 0.07.

Although somewhat more dispersed, the data for the square tube also follow a power law with an exponent -0.67 ± 0.10 compatible with -2/3. According to Eq. (3.1), we rescale R by the square side c, and we plot R/c as a function of v in the same plot as the data for circular tubes (Fig. 5). All data collapse on the same master curve, showing that the hydraulic diameter is indeed the relevant geometrical length to compare tubes of different cross-sectional shapes. This master curve is a power law of exponent -0.69±0.03, again compatible with an exponent -2/3.

Concerning the prefactor, the best fit of the data with a power law of exponent -2/3

gives Rv 2/3 /b = 0.74 ± 0.02 m 2/3 , whereas inserting the measured solution properties in the prediction (3.6) gives Rv 2/3 /b = 2.08 m 2/3 for β = 1 (mobile interfaces), and 1.66 m 2/3 for β = 2 1/3 (rigid interfaces).

Unsteady behaviour and rupture of the soap films

We now test the predictions of the minimal radius of curvature given by (3.2) and (3.4).

The agreement with the most curved films that we were able to detect is excellent, within 2% for the square and 6% for the circular tubes (Fig. 5). No steady films have been recorded at higher velocity: they broke before entering the camera field of view.

To go beyond the limit of stability, we produced films at low velocity and accelerated them afterwards, typically 10 cm before they enter the field of view. With this procedure, when the gas flux is suddenly increased, the film is accelerated faster than the PBs, thus becomes much drier. Many films burst during this transient stage (Fig. 2c), but not all. With this procedure, few points have been obtained above the theoretical velocity limit on Fig. 5. They are still compatible with the friction law curve and their radius of curvature are then smaller than the theoretical predictions. As depicted on Fig. 2d, a long cylindrical neck is then present at the rear of the film. In order to evaluate the stability of these surprising shapes, we recorded images with a larger field of view (15 cm). These shape appeared to be unsteady. Several scenarii have been observed: breakage after few centimeters; evolution towards the half sphere when the velocity is very close to the critical velocity; pinch-off of the neck into a separate bubble in some cases, that does not touch the wall anymore.

In contrast, we also checked that the shape recorded at velocity below the critical velocity does not change during the 15 cm long trajectory.

Discussion

The weak influence of the Plateau border size

We have hitherto disregarded the role of the PB size, which does not appear in Eq. (3.6) with which our measurements have been interpreted. However, we have observed in many cases (e.g. Fig. 2b) that the soap films are not symmetric. The top part of the PB, which is also the driest because of drainage, lags often significantly behind the bottom part, hence experiences more friction.

To quantify this possible influence of the PB size, we have extracted for the experiments with the circular tube of diameter 5.9 mm the two angles, ϕ top and ϕ bottom , with which the soap film meets the top and bottom part of the tube. This gives a ratio on the friction force per unit PB length: f (v) top /f (v) bottom = sin ϕ top / sin ϕ bottom . For all experiments but one, we have checked that log |f (v) top /f (v) bottom | remains below 0.1, which is the order of magnitude of the scatter of the data points with respect to the master curve in Fig. 5.

Therefore, the PB size is a correction to the main variation of the radius of curvature, with the velocity, that can be neglected in a first approach. In the absence of a theoretical prediction, it may be possible to include its role using empirical laws such as the one proposed by [START_REF] Raufaste | Dissipation in quasi-two-dimensional flowing foams[END_REF].

Boundary condition: rigid or mobile interfaces?

SDS is known to be very mobile at an air/water interface and to exchange quickly with the solution bulk in case of area variation. Nevertheless, the expansion rate in our experiment scales like v/ℓ ∼ v/rCa 1/3 , with ℓ the length of the transition domain, and reaches 10 5 s -1

for the highest velocities. At such a high rate, concentration gradients can appear, which effect is to impede surface area variations, through Marangoni stress. The boundary condition to be used is thus not obvious. A local conservation of the interface area, as used by [START_REF] Denkov | Wall slip and viscous dissipation in sheared foams: Effect of surface mobility[END_REF] for finite size bubbles, imposes that the whole wetting film interface moves at the film velocity v and is thus unacceptable : in our geometry, the total interface area decreases downstream and increases upstream. In Appendix A, we allow the wetting film interface to move with the wall, as a rigid interface. The area variation is assumed to occur on larger length scale and smaller time scale in the moving film or in the PB interfaces and to be without influence on f (v). The viscous force for a rigid interface is 1.26 times larger than for a mobile interface, but it is still lower than our experimental result.

Contribution of the Plateau border to the viscous force

In Bretherton's classical result, the friction force is localised in the transition zone between the wetting film and an external infinite reservoir. Here, on the contrary, the PB has a finite size r. Recirculations, with velocities larger than the film velocity have been observed in the PB, when tiny bubbles were trapped there; they have also been predicted in numerical simulations [START_REF] Saugey | Wall slip of bubbles in foams[END_REF]). This induces extra velocity gradients, estimated in Appendix B, showing that the contribution of the PB to the viscous force should be of the same order than the contribution of the transition domains. We show that this additional contribution scales as the Bretherton term, which could explain why the scaling f (v) ∼ σCa 2/3 still holds in our experiments.

Figure 1 .

 1 Figure 1. Sketch of the experimental setup.

Figure 2 .

 2 Figure 2. Snapshots of soap films pushed from left to right through a tube of outer diameter 9.0 mm, at velocities (a) 0.37, (b) 0.64, (c) 1.81 and (d) 2.47 m/s. Films (a) and (b) were observed to be stable, and films (c) and (d) displayed an unsteady behaviour: film (c) broke down as shown in the three consecutive displayed frames (interframe 0.33 ms), whereas film (d) tends to detach from its Plateau border (PB) through a long cylindrical neck, which eventually leads to film bursting or pinch-off of a bubble (see Sec. 4.2). Notice the top-bottom asymmetry of the PB, due to drainage [the direction of gravity is indicated in (a)], and, in (a) and (b), the damped oscillations of the wetting film ahead of the PB in contrast with the monotonous profile of the wetting film at the rear of the PB.

Figure 3 .

 3 Figure 3. Snapshots of soap films pushed from left to right through a square tube of side 10.0 mm, (a) in tilted view, to get an impression of the whole shape, and (b) facing one of the four sides: the front part of the shape is the trace of the central part of the film, and the back part is the PB in contact with the imaged side. This is the view used to measure the radius of curvature of the film. The superimposed notations are used in Sec. 3.

Figure 4 .

 4 Figure 4. (a) Example of a spatiotemporal diagram. The film is represented by a straight dark band. To measure the film curvature, the raw image (b) is binarised (c), and (d) the coordinates of the film contour are extracted (dotted curve) and fitted by a circle (dashed curve), also superimposed in (b) with various notations used in the text.

  .2) and f (v) cannot exceed 2σ.For a square tube, the film shape should be of constant curvature and compatible with the local viscous friction on the boundary. Let z = ζ(x, y) be the equation of the film in the square domain 0 ≤ x, y ≤ a. We consider an elementary portion of a PB in contact B.Dollet and I. Cantat

  where n = (-ζ x e x -ζ y e y + e z )/ 1 + ζ 2 x + ζ 2 y is the normal unit vector of the transverse film, and t ′ = (e x + ζ x e z )/ 1 + ζ 2 x is the tangential unit vector along the portion of PB. Hence, t • e z = ζ y / (1 + ζ

  to predict the viscous dissipation associated to the flow of foams against a wall. Finally, combining (3.1) and (3.5), we get the prediction for the radius of curvature of the film:

Figure 5 .

 5 Figure 5. Plot of log R/b as a function of log v (expressed in m/s), for circular tubes of inner diameters 3.6 (•), 5.9 ( ) and 8.8 mm (♦), and for the square tube of side 10.0 mm (×). The three points with the symbol were obtained in unsteady situations, see Sec. 4.2. The plain line represents the best linear fit through all data; its equation is log R/a = (-0.69 ± 0.03) log v -0.131 ± 0.008. The dashed (dotted) line represents the theoretical limit of stability for soap films in circular (square) tubes. Inset: plot of log R (in mm)as a function of log v, and its best linear fit, for the tube of diameter 8.8 mm.

Absence of inertial correction

Inertial effects can be estimated using the Weber number W = ρv 2 r/σ, with r ∼ 10 -4 m the PB characteristic size and ρ = 10 3 kg/m 3 the solution density. For our velocity range, W varies between 0.04 and 25: inertial effects should then be visible and even dominant in that velocity range. The viscous force scaling in Ca 2/3 , signature of the visco-capillary regime, is therefore somewhat unexpected in this regime. Actually, the friction force scales as f (v) ∼ ηℓv/h ∞ . [START_REF] Quéré | Le mouillage dynamique des fibres[END_REF] have considered the deviations of h ∞ and ℓ from the visco-capillary scaling as inertia becomes significant, see Eq. (3.20) in their article. However, the scaling ℓ/h ∞ ∼ Ca -1/3 remains whatever the value of W .

Hence, since inertia modifies ℓ and h ∞ likewise, the scaling f (v) ∼ σCa 2/3 is unaffected by the inertial corrections, as observed experimentally here.

Concluding remarks

Our results show that the limit of stability of soap films pushed through tubes obeys simple theoretical predictions, which match the experiments with an excellent precision.

Once the friction law is known, this gives a quantitative criterion on the maximal velocity of a foam flow in a confined geometry.

Concerning the viscous friction between the soap film and the wall, the scaling proposed by [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF] was shown to hold beyond its expected range of validity, at rather high capillary and Weber numbers. In the frame of the long-standing debate on the boundary condition and the role of surfactants [START_REF] Ratulowski | Marangoni effects of trace impurities on the motion of long gas bubbles in capillaries[END_REF][START_REF] Shen | Fiber coating with surfactant solutions[END_REF][START_REF] Denkov | Wall slip and viscous dissipation in sheared foams: Effect of surface mobility[END_REF], this calls for further investigations on the viscous dissipation associated with the motion of bubbles, films and foams close to walls. In this context, we are currently studying the role of surface viscoelasticity on soap film and foam flows in confined geometries.

Appendix A. Friction law: mobile versus rigid interfaces

Here we compute the contribution to f of the two transition domains between the PB and the wetting film, upstream (top sign in the equations when the symbols ± or ∓ appear) and downstream (bottom sign in the equations). We use the frame of the steady film in which the wall moves at velocity v. Both contributions f ± are determined independently and in both cases we orient the x axis towards the thin wetting film. The velocity v is thus positive for the upstream interface and negative for the downstream interface. The lubrication equation allows to determine the velocity field as a function of the interface profile h(x): u(x, y) = (σ/η) λhy -1 2 y 2 ∂ 3 h/∂x 3 + v, with λ = 1 if the tangential stress is vanishing at the interface and λ = 1/2 if the interface moves at the wall velocity v. The governing equation for the interface profile is, with Ca = η|v|/σ and β = 2/(3λ -1), ±3βCa(h ∞ -h) = h 3 ∂ 3 h/∂x 3 . Using the rescaling h = Hh ∞ and x = Xh ∞ /(3βCa) 1/3 with h ∞ the wetting film thickness we get H ′′′ H 3 = ±(1-H), with boundary conditions H(∞) = 1, H ′ (∞) = 0 and H ′′ (∞) = 0, where the prime denotes derivation with respect to X. Finally, the contribution to the viscous force is given by

∂u/∂y(y = 0)dx that can be expressed as a function of the unknown solution H ± (X) and of the various parameters involved in the rescaling. We obtain:

Finally the contribution f ± of the two transition domains verifies

26. This ratio is twice smaller than the thickening factor h ∞ (rigid)/h ∞ (mobile) = 4 2/3 shown by [START_REF] Ratulowski | Marangoni effects of trace impurities on the motion of long gas bubbles in capillaries[END_REF] and [START_REF] Shen | Fiber coating with surfactant solutions[END_REF]. Actually, the wetting film thickness is determined Hence, h ∞ and f scale differently on λ and β, whence their different ratios between the mobile and rigid limits.

Appendix B. Contribution of the Plateau border to the friction

In addition to the contribution f ± of the transition domains, the PB itself has a nonnegligible contribution f P B to the total viscous force per unit length. For a film flowing at velocity v in a tube of radius a, with a PB of typical size r, this contribution scales like f P B ∼ ηrv P B /r ∼ ηv P B .

v P B can either be governed by the wall motion, which implies v P B ∼ v, or by a pressure gradient due to the fore-aft asymmetry of the PB (Fig. 6), which implies v P B ∼ (P + -P -)r/η. We compare both terms hereafter. The pressures are determined from the geometrical constraints governing the PB shape (see Fig. 6). Using the Laplace laws P ± g -P ± = σ/r ± , P + g -P - g = 4σ sin φ/a and the geometrical relation r ± (1∓ sin φ) = r, we get P + -P -∼ (P + g -P - g )(1 + a/2r) ∼ a(P + g -P - g )/2r. With f = a(P + g -P - g )/2

we get finally f P B /f ∼ 1. The dissipation in the PB is thus not negligible in comparison with that in the transition region, and has the same scaling. From f P B ∼ ηv P B we deduce Ca P B ∼ Ca 2/3 , or v P B /v ∼ Ca -1/3 > 1. The velocity in the PB is thus mainly governed by pressure effects and not only by the boundary velocity.


Pascal Koiran

Natacha Portier

Sébastien Tavenas
email: .tavenas@ens-lyon.fr.

Stéphan Thomassé Lip

A τ -Conjecture for Newton Polygons

One can associate to any bivariate polynomial P (X, Y ) its Newton polygon. This is the convex hull of the points (i, j) such that the monomial X i Y j appears in P with a nonzero coefficient. We conjecture that when P is expressed as a sum of products of sparse polynomials, the number of edges of its Newton polygon is polynomially bounded in the size of such an expression. We show that this "τ -conjecture for Newton polygons," even in a weak form, implies that the permanent polynomial is not computable by polynomial size arithmetic circuits. We make the same observation for a weak version of an earlier "real τ -conjecture." Finally, we make some progress toward the τ -conjecture for Newton polygons using recent results from combinatorial geometry.

Introduction

Let f ∈ Z[X] be a univariate polynomial computed by an arithmetic circuit (or equivalently, a straight-line program) of size s starting from the variable X and the constant 1. According to Shub and Smale's τ -conjecture [START_REF] Shub | On the intractability of Hilbert's Nullstellensatz and an algebraic version of "P=NP[END_REF], the number of integer roots of f should be bounded by a fixed polynomial function of s. It was shown in [START_REF] Shub | On the intractability of Hilbert's Nullstellensatz and an algebraic version of "P=NP[END_REF] that the τ -conjecture implies a P = NP result for the Blum-Shub-Smale model of computation over the complex numbers [START_REF] Blum | On a theory of computation and complexity over the real numbers: NP-completeness, recursive functions and universal machines[END_REF][START_REF] Blum | Algebraic settings for the problem "P =NP?[END_REF]. A similar result was obtained by Bürgisser [START_REF] Bürgisser | On defining integers and proving arithmetic circuit lower bounds[END_REF] for another algebraic version of P versus NP put forward by Valiant [START_REF] Valiant | Completeness classes in algebra[END_REF][START_REF] Valiant | Reducibility by algebraic projections[END_REF] at the end of the 1970's. A succinct way of stating this VP versus VNP problem goes as follows: can we compute the permanent of a n × n matrix with a number of arithmetic operations which is polynomial in n? This question can be formalized using the computation model of arithmetic circuits. The permanent plays a special role here because it is VNP-complete, and it can be replaced by any other VNP-complete family of polynomials. We refer to Bürgisser's book [START_REF] Bürgisser | Completeness and Reduction in Algebraic Complexity Theory[END_REF] for an introduction to this topic and to two recent surveys on arithmetic circuit complexity by Shpilka and Yehudayoff [START_REF] Shpilka | Arithmetic circuits: A survey of recent results and open questions[END_REF] and by Chen, Kayal and Wigderson [START_REF] Chen | Partial derivatives in arithmetic complexity and beyond[END_REF].

As a natural approach to the τ -conjecture, one can try to bound the number of real roots instead of the number of integer roots. This fails miserably since the number of real roots of a univariate polynomial can be exponential in its arithmetic circuit size. Chebyshev polynomials provide such an example [START_REF] Smale | Mathematical problems for the next century[END_REF] (a similar example was provided earlier by Borodin and Cook [START_REF] Borodin | On the number additions to compute specific polynomials[END_REF]). A real version of the τ -conjecture was nevertheless proposed in [START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF]. In order to avoid the aforementioned counterexamples, the attention is restricted to arithmetic circuits of a special form: the sums of products of sparse polynomials. In spite of this restriction, the real τ -conjecture still implies that the permanent is hard to compute for general arithmetic circuits [START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF].

In this paper, we propose a τ -conjecture for Newton polygons of bivariate polynomials. Like the real τ -conjecture, it deals with sums of products of sparse polynomials and implies that the permanent is hard for general arithmetic circuits. A common idea to these three τ -conjecture is that "simple" arithmetic circuits should compute only "simple" polynomials. In the original τ -conjecture, the simplicity of a polynomial is measured by the number of its integer roots; in the real τ -conjecture it is measured by the number of its real roots; and in our new conjecture by the number of edges of its Newton polygon. These conjectures are independent in the sense that we know of no implication between two of them.

Organization of the paper

In Section 2 we review some basic facts about Newton polygons and formulate the corresponding τ -conjecture. We also state in Theorem 1 the motivating result for this paper: a proof of the conjecture, even in a very weak form, implies a lower bound for the permanent. In Section 3 we give a proof of this result and of a refinement: it suffices to work with sums of powers of sparse polynomials rather than with sums of arbitrary products. We also point out that this refinement applies to the real τ -conjecture from [START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF], and that (like in Theorem 1) a very weak form of this conjecture implies a lower bound for the permanent. These observations improve the results stated in [START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF]. In Section 4 we use a recent result of convex geometry [START_REF] Eisenbrand | Convexly independent subsets of the Minkowski sum of planar point sets[END_REF] to provide nontrivial upper bounds on the number of edges of Newton polygons. Our results fall short of establising the new τ -conjecture (even in the weak form required by Theorem 1) but they improve significantly on the naive bound obtained by brute-force expansion. For instance, as a very special case of our results we have that the Newton polygon of f g + 1 has O(t 4/3 ) edges if the bivariate polynomials f and g have at most t monomials. The straightforward bound obtained by expanding the product f g is only O(t 2 ). We conclude the paper with a couple of open problems. In particular, we ask whether this O(t 4/3 ) upper bound is optimal. In the appendix, we improve on this upper bound by giving a linear upper bound in a special case.

Newton Polygons

In the following we will consider a field K of characteristic 0. In fact, most of the results of this paper are valid for a general field.

We first recall some standard background on Newton polygons. Much more can be found in the survey [START_REF] Sturmfels | Polynomial equations and convex polytopes[END_REF]. Consider a bivariate polynomial f ∈ K[X, Y ]. To each monomial X i Y j appearing in f with a nonzero coefficient we associate the point with coordinates (i, j) in the Euclidean plane. We denote by Mon(f ) this finite set of points. If A is a set of points in the plane, we denote by conv(A) the convex hull of A. By definition the Newton polygon of f , denoted Newt(f ), is the convex hull of Mon(f ) (In particular, Newt(f ) = conv(Mon(f ))). Note that Newt(f ) has at most t edges if f has t monomials. It is well known [START_REF] Ostrowski | ber die Bedeutung der Theorie der konvexen Polyeder fr die formale Algebra[END_REF] that the Newton polygon of a product of polynomials is the Minkowski sum of their Newton polygons, i.e., Newt(f g) = Newt(f ) + Newt(g) = {p + q; p ∈ Newt(f ), q ∈ Newt(g)}.

A short proof of this fact can be found in [START_REF] Gao | Absolute irreducibility of polynomials via Newton polytopes[END_REF]. As a result, if f has s monomials and g has t monomials then Newt(f g) has at most s + t edges. More generally, for a product f = g 1 g 2 • • • g m , Newt(f ) has at most m i=1 t i edges where t i is the number of monomials of g i ; but f can of course have up to m i=1 t i monomials. The number of edges of a Newton polygon is therefore easy to control for a product of polynomials. In contrast, the situation is not at all clear for a sum of products. We propose the following conjecture.

Conjecture 1 (τ -conjecture for Newton polygons).

There is a polynomial p such that the following property holds.

Consider any bivariate polynomial

f ∈ K[X, Y ] of the form f (X, Y ) = k i=1 m j=1 f ij (X, Y ) (1)
where the f ij have at most t monomials. Then the Newton polygon of f has at most p(kmt) edges.

Note that the bound only depends on k, m and t, and so does not depend on the degrees and coefficients of the polynomials f i,j .

The "real τ -conjecture" [START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF] is a similar conjecture for real roots of sums of products of sparse1 univariate polynomials, and it implies that the permanent does not have polynomial-size arithmetic circuits. We are not able to compare these two conjectures. However, as we shall see shortly, Conjecture 1 also implies that the permanent does not have polynomial-size arithmetic circuits.

One can also formulate a τ -conjecture for the multiplicities of (complex) nonzero roots of univariate polynomials (see Corollary 2.4.1 in [START_REF] Hrubeš | On the Real τ -Conjecture and the Distribution of Complex Roots[END_REF]). This conjecture again implies that the permanent is hard for arithmetic circuits.

By expanding the products in [START_REF] Agrawal | Arithmetic circuits: A chasm at depth four[END_REF] we see that f has at most k.t m monomials, and this is an upper bound on the number of edges of its Newton polygon. In order to improve this very coarse bound, the main difficulty we have to face is that the k-fold summation in the definition of f may lead to cancellations of monomials. As a result, some of the vertices of Newt(f ) might not be vertices of the Newton polygons of any of the k products m j=1 f ij (X, Y ). We give two examples of such cancellations below. If there are no cancellations (for instance, if the f ij only have positive coefficients) then we indeed have a polynomial upper bound. In this case, Newt(f ) is the convex hull of the union of the Newton polygons of the k products. Each of these k Newton polygons has at most mt vertices, so Newt(f ) has at most kmt vertices and as many edges.

Example 1. We define A(X, Y ) = XY +X 2 +X 2 Y 2 +X 3 Y +X 5 Y , B(X, Y ) = 1 + XY 2 , C(X, Y ) = -X -XY -X 2 Y 2 and D(X, Y ) = Y + X + X 2 Y + X 4 Y . AB point of Mon(AB) CD point of Mon(CD) Then, AB + CD =(XY + X 2 + X 2 Y 2 + X 3 Y + X 5 Y + X 2 Y 3 + X 3 Y 2 + X 3 Y 4 + X 4 Y 3 + X 6 Y 3 ) -(XY + X 2 + X 3 Y + X 5 Y + XY 2 + X 2 Y + 2X 3 Y 2 + X 5 Y 2 + X 2 Y 3 + X 4 Y 3 + X 6 Y 3 ) =X 2 Y 2 + X 3 Y 4 -XY 2 -X 2 Y -X 3 Y 2 -X 5 Y 2 AB + CD AB AB + CD CD point of Mon(AB + CD)
The two rectangle points lie on the convex hull of Mon(AB + CD), but do not lie on the convex hulls of Mon(AB) or Mon(CD).

Example 2. We define f (X, Y ) = 1 + X 2 Y + Y 2 X, g(X, Y ) = 1 + X 4 Y + XY 4
and we consider Mon(f g -1).

point of Mon(f g)

The two rectangle points lie on the convex hull of Mon(f g -1), but do not lie on the convex hull of Mon(f g).

Conjecture 1 implies that the permanent is hard for arithmetic circuits. In fact, a significantly weaker bound on the number of edges would be sufficient: Theorem 1. Assume that for some universal constant c < 2, the upper bound 2 (m+log kt) c on the number of edges of Newt(f ) holds true for polynomials of the form (1) whenever the product kmt is sufficiently large. Then the permanent is not computable by polynomial size arithmetic circuits.

For instance, an upper bound of the form 2 O(m) (kt) O(1) would be sufficient. Note that the parameter m plays a very different role than the parameters k and t.

Remark 1. Theorem 1 holds true even in the case of a field of positive characteristic different from two. In characteristic two, the permanent is equal to the determinant and is therefore easy for arithmetic circuits. However, the conclusion of Theorem 1 remains true in characteristic two if we replace the permanent by the Hamiltonian polynomial (or any other polynomial family which is VNP-complete in characteristic two).

Proof of Theorem 1, and a Refinement

Consider the polynomial

f n (X, Y ) = 2 n -1 i=0 X i Y i 2 .
(

) 2
The Newton polygon of f n has exactly 2 n edges. Our proof of Theorem 1 is by contradiction. Assuming that the permanent is computable by polynomial size arithmetic circuits we will show that f n can be put under form [START_REF] Agrawal | Arithmetic circuits: A chasm at depth four[END_REF]

with k = n O( √ n log n) , t = n O( √ n log n) and m = O( √ n).
Note that the upper bound on m is much smaller than those on k and t. Then, from the assumption in Theorem 1 we conclude that Newt(f n ) has at most 2 (m+log kt) c edges. This is a contradiction since for large enough n, this upper bound is smaller than the actual number of edges of Newt(f n ), namely, 2 n (here, we use the fact that the constant c in Theorem 1 is smaller than 2).

Reduction of arithmetic circuits to depth 4 is an important ingredient in the proof of Theorem 1. This phenomenon was discovered by Agrawal and Vinay [START_REF] Agrawal | Arithmetic circuits: A chasm at depth four[END_REF]. We will use it under the following form [START_REF] Koiran | Arithmetic circuits: the chasm at depth four gets wider[END_REF] (recall that a depth 4 circuit is a sum of products of sums of products of inputs; sum and product gates may have arbitrary fan-in). Note that Theorem 3 of [START_REF] Koiran | Arithmetic circuits: the chasm at depth four gets wider[END_REF] provides this bound for the case where d is the so-called "formal degree" of C rather than the degree of the polynomial computed by C. Theorem 2 as stated above can then be derived by an application of the standard homogenization trick (see e.g. Proposition 5 and Theorem 5 in [START_REF] Koiran | Arithmetic circuits: the chasm at depth four gets wider[END_REF]). It was recently shown [START_REF] Tavenas | Improved bounds for reduction to depth 4 and depth 3[END_REF] that the size bound for Γ can be reduced from t O(

√ d log d) to t O( √ d)
when d is polynomially bounded in t; this improvement preserves the O( √ d) bound on the fan-in of multiplication gates. We can now complete the proof of Theorem 1. The idea of this proof is based on the one of Proposition 2 in [START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF], but is simpler since it doesn't use the counting hierarchy. A similar simplification can be done for the real τconjecture [START_REF] Tavenas | Bornes inférieures et supérieures pour les circuits arithmétiques[END_REF].

We can expand the exponents i and i 2 of (2) in base 2. This leads to the equality

f n (X, Y ) = h n (X 2 0 , X 2 1 , . . . , X 2 n-1 , Y 2 0 , . . . , Y 2 2n-1 ) (3)
where h n (x 0 , x 1 , . . . , x n-1 , y 0 , . . . , y 2n-1 ) is the multilinear polynomial

α∈{0,1} n ,β∈{0,1} 2n a(n, α, β)x α 0 0 x α 1 1 • • • x α n-1 n-1 y β 0 0 • • • y β 2n-1 2n-1 .
Here the exponents α j , β j denote the coordinates of the vectors α, β. The coefficient a(n, α, β) equals 1 if the vector β corresponds to the square of the vector α and 0 otherwise. Note that h n is a polynomial in 3n variables. The coefficients a(n, α, β) can be computed in time polynomial in n. By Valiant's criterion [START_REF] Bürgisser | Completeness and Reduction in Algebraic Complexity Theory[END_REF], this implies that the polynomial family (h n ) belongs to the complexity class VNP. Since the permanent is VNP-complete and is assumed to have polynomial-size circuits, (h n ) also has polynomial-size circuits. By Theorem 2, it follows that the polynomials h n are computable by depth 4 circuits of size n O( √ n log n) with multiplication gates of fan-in O( √ n). Using (3), we can plug in powers of X, Y and powers of 2 into those circuits to express f n as a sum of products like in [START_REF] Agrawal | Arithmetic circuits: A chasm at depth four[END_REF]. The resulting parameters k and t are of order n O( √ n log n) , and m = O( √ n). As explained at the beginning of this section, this leads to a contradiction with the assumption in Theorem 1.

In the remainder of this section we give a refinement of Theorem 1. We show that it suffices to bound the number of edges of the Newton polygons of sums of powers of sparse polynomials in order to obtain a lower bound for the permanent.

Theorem 3. Fix a universal constant c < 2, and assume that we have the upper bound 2 (m+log kt) c on the number of edges of Newt(f ) for polynomials of the form

f (X, Y ) = k i=1 a i f i (X, Y ) m (4)
where a i ∈ K and the f i have at most t monomials (as in Theorem 1, we require this upper bound to hold only if kmt is sufficiently large). Then the permanent is not computable by polynomial size arithmetic circuits.

Clearly, we can assume that all the coefficients a i are equal to 1 (multiply f i by a m-th root of a i if necessary).

Theorem 3 is an easy consequence of Theorem 1 and Fischer's formula [START_REF] Fischer | Sums of like powers of multivariate linear forms[END_REF]. This formula shows that any monomial z 1 z 2 • • • z m can be expressed as a linear combination of 2 m-1 powers of linear forms.

Lemma 1. For any m, we have

2 m-1 m!z 1 z 2 • • • z m = r=(r 2 ,...,rm)∈{-1,1} m-1 ( m i=2 r i )(z 1 + m i=2 r i z i ) m .
Note that the exponential blowup entailed by Fischer's formula is acceptable because we will apply it with a value of m which is small compared to the main complexity parameter n, i.e., with m = O( √ n). The idea of using Fischer's formula to turn a product into sums of powers comes from [START_REF] Gupta | Arithmetic circuits: A chasm at depth three[END_REF][START_REF]An exponential lower bound for the sum of powers of bounded degree polynomials[END_REF].

Proof of Theorem 3. We show that the assumption in Theorem 3 implies that of Theorem 1. Consider therefore a polynomial f of the form (1). We rewrite it as a sum of powers by applying Lemma 1 to each of the k products in (1). This yields an identity of the form

f (X, Y ) = k ′ i=1 a i f i (X, Y ) m
where a i ∈ K, the f i have at most mt monomials, and k ′ = 2 m-1 k. We are now in position to apply Theorem 1: Newt(f ) has at most 2 (m+log k ′ t ′ ) c edges. For any constant c ′ > c, this is less than 2 (m+log kt) c ′ if kmt is sufficiently large. We have therefore derived the hypothesis of Theorem 1 from that of Theorem 3, and we can conclude that the permanent is hard for arithmetic circuits.

Remark 2. For Theorem 3, the fact that K has characteristic 0 is important. Indeed, in positive characteristic Lemma 1 does not allow to rewrite a monomial z 1 z 2 • • • z m as a linear combination of powers of linear forms. Remark 3. As pointed out in the introduction, we gave in [START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF] similar results for real roots of univariate polynomials rather than for Newton polygons of bivariate polynomials. More precisely, let us measure the size of a sum of products of sparse polynomials by s = kmt. This definition of "size" applies to bivariate polynomials of the form (1) as well as to their univariate analogues. We showed that for any constant c < 2, a 2 (log s) c upper bound on the number of real roots implies that the permanent is hard for arithmetic circuit (see Conjecture 3 in [START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF] and the remarks following it). In fact, the same proof shows than an upper bound of the form 2 (m+log kt) c as in Theorem 1 is sufficient. This is clearly a better way of stating our result since it allows for a much worse dependency of the number of real roots with respect to m. Moreover, as in Theorem 3 it is sufficient to establish this bound for sums of powers. As in the proof of Theorem 3, this follows from a straightforward application of Fischer's formula.

Upper Bounds from Convexity Arguments

In this section we improve the coarse upper bound k.t m upper bound on the number of edges of Newton polygons of polynomials of the form (1). These results apply to fields of arbitrary characteristic. Our main tool is a result of convex geometry [START_REF] Eisenbrand | Convexly independent subsets of the Minkowski sum of planar point sets[END_REF]. Theorem 4. Let P and Q be two planar point sets with |P | = s and |Q| = t. Let S be a subset of the Minkowski sum P +Q. If S is convexly independent (i.e., its elements are the vertices of a convex polygon) we have |S| = O(s 2/3 t 2/3 + s + t).

It is known that this upper bound is optimal up to constant factors [START_REF] Bílka | A tight lower bound for convexly independent subsets of the Minkowski sums of planar point sets[END_REF] (a non-optimal lower bound was also given in [START_REF] Swanepoel | Large convexly independent subsets of Minkowski sums[END_REF]). A linear upper bound is known [START_REF] Halman | The convex dimension of a graph[END_REF] for the case where P and Q are convexly independent.

We first consider sums of products of two polynomials.

Theorem 5. Consider a bivariate polynomial f ∈ K[X, Y ] of the form

f (X, Y ) = k i=1 f i g i (X, Y ) (5)
where the f i have at most r monomials and the g i have at most s monomials. The Newton polygon of f has O(k(r 2/3 s 2/3 + r + s)) edges.

Proof. Let S i be the set of points in the plane corresponding to the monomials of f i g i which appear in f with a nonzero coefficient. Since Newt(f ) is the convex hull of k i=1 conv(S i ), it is enough to bound the number of vertices of conv(S i ). Those vertices form a convexly independent subset of the Minkowski sum Mon(f i )+Mon(g i ). By Theorem 4, it follows that conv(S i ) has O(r 2/3 s 2/3 + r + s) vertices. Multiplying this estimate by k yields an upper bound on the number of vertices and edges of Newt(f ).

From this result it is easy to derive an upper bound for the general case, where we have products of m ≥ 2 polynomials. We just divide the m factors into two groups of approximately m/2 factors, and in each group we expand the product by brute force.

Theorem 6. Consider any bivariate polynomial

f ∈ K[X, Y ] of the form f (X, Y ) = k i=1 m j=1 f ij (X, Y ) (6)
where m ≥ 2 and the f ij have at most t monomials. The Newton polygon of f has O(k.t 2m/3 ) edges.

Proof. As suggested above, we write each of the k products as a product of two polynomials

F i = ⌊m/2⌋ i=1 f i and G i = ⌈m/2⌉ i=1 f i .
We can now apply Theorem 5 to f = k i=1 F i G i , with r = t ⌊m/2⌋ and s = t ⌈m/2⌉ . In the resulting O(k(r 2/3 s 2/3 + r + s)) upper bound the term kr 2/3 s 2/3 dominates since r 2/3 s 2/3 = t 2(⌊m/2⌋+⌈m/2⌉)/3 = t 2m/3 and m ≥ 2.

In order to avoid the brute force expansion in the proof of this theorem it is natural to consider for each i a convexly independent subset S i of the Minkowski sum of the m sets Mon(f i1 ), . . . , Mon(f im ). This is exactly the open problem at the end of [START_REF] Bílka | A tight lower bound for convexly independent subsets of the Minkowski sums of planar point sets[END_REF]: determine the maximal cardinality M m (t) of a convexly independent subset of the Minkowski sum of m sets P 0 , . . . , P m-1 of t points in the Euclidean plane. For instance, the lower bound of [START_REF] Bílka | A tight lower bound for convexly independent subsets of the Minkowski sums of planar point sets[END_REF] combined with the upper bound of [START_REF] Eisenbrand | Convexly independent subsets of the Minkowski sum of planar point sets[END_REF] shows that M 2 (t) = Θ(t 4/3 ). Unfortunately, we shall see that M m (t) = t Ω(m) , so that brute force expansion is not very far from the optimum. More explicitly,

P k = {(b 2k .i, b k .j); 0 ≤ i ≤ b 2 -1 and 0 ≤ j ≤ b -1}.
Clearly, the Minkowksi sum P 0 + . . .

+ P m-1 is the grid {0, . . . , b 2m -1} × {0, . . . , b m -1} of size b 2m × b m .
The next lemma (which is certainly not optimal) shows how to find a fairly large set of convexly independent points in a grid. Lemma 2. If n(n -1)/2 < M and n < N it is possible to find n convexly independent points in the grid M × N .

Proof. We start from the origin and build a sequence of n -1 line segments. The i-th segment has horizontal length i and slope 1/i. We can keep going as long as we do not go out of the grid, i.e., as long as n(n -1)/2 < M and n < N . Altogether, the n -1 segments have n endpoints and they are convexly independent.

Proposition 1. For all m and infinitely many values of t we have: This result shows that other ingredients than Theorem 4 will be needed to answer Conjecture 1 positively. A similar argument can be made for the case where the sets P 0 , P 1 , . . . , P m-1 in the Minkowski sum are all equal (this is a natural case to look at in light of Theorem 3, which shows that it suffices to deal with sums of powers in order to obtain a lower bound for the permanent). More precisely, let M ′ m (t) be the maximal cardinality of a convexly independent subset of an m-fold Minkoski sum P + P + • • • + P where P is a set of at most m points. By definition we have M ′ m (t) ≤ M m (t). In the other direction we have M ′ m (t) ≥ M m (⌊t/m⌋): just replace the m sets of size ⌊t/m⌋ by their union. Hence we have M ′ m (t) ≥ ⌊t/m⌋ m/3 -1.

M m (t) ≥ t m/3 -1.

Final Remarks

In this paper we have proposed a conjecture on the number of edges of the Newton polygon of a sum of products of sparse polynomials; and we have shown in Theorem 1 that even a weak version of this conjecture implies a lower bound for the permanent. We conclude with a couple of additional open problems.

1. Consider two polynomials f, g ∈ K[X, Y ] with at most t monomials each. What is the maximum number of edges on the Newton polygon of f g + 1?

The difficult case appears when the constant term of f g equals -1.

Theorem 5 provides a O(t 4/3 ) upper bound, but as far as we know the "true" bound could be linear in t. In the appendix we prove a linear upper bound under the assumption that f and g have the same supports (i.e., Mon(f ) = Mon(g)) and that the square of any nonconstant monomial appearing in f and g does not appear in f g.

2.

More generally, what is the maximum number of edges on the Newton polygon of f 1 . . . f m + 1, where the f i again have at most t monomials? Theorem 6 provides a O(t 2m/3 ) upper bound, but the true bound could be of the form 2 O(m) t O (1) ; it could even be polynomial in m and t, as implied by Conjecture 1.

Theorem 2 .

 2 Let C be an arithmetic circuit of size t computing a polynomial of degree d. There is an equivalent depth four circuit Γ of size t O( √ d log d) with multiplication gates of fan-in O( √ d).

Example 3 .

 3 Fix an integer b ≥ 2. Let P k be the b 2 × b grid made of the integer points (x, y) such that: -all the digits in base b 2 of x are equal to zero, except possibly the digit of weight b 2k ; -all the digits in base b of y are equal to zero, except possibly the digit of weight b k .

Proof.

  From Example 3 and Lemma 2 (choosing M = b 2m and N = b m ) we have M m (b 3 ) ≥ n if n(n -1) < 2b 2m and n < b m . Hence M m (b 3 ) ≥ b m -1. The result follows by setting t = b 3 .

Here and in[START_REF] Koiran | Shallow circuits with high-powered inputs[END_REF], the term "sparse" refers to the fact that we measure the size of a polynomial fij by the number of its monomials.

Acknowledgments

Proposition 1 arose from a discussion with Mark Braverman and an improvement was suggested by an anonymous referee. We thank the three referees for suggesting several improvements in the presentation of the paper.

Appendix: the Newton polygon of f g + 1

In this section we denote by 0 the point in the plane with coordinates (0, 0).

We give here (in Theorem 7) a linear upper bound assuming the following two properties:

(i) The polynomials f and g have the same support, i.e., Mon(f ) = Mon(g). We denote by {p 0 , . . . , p t-1 } this common support.

(ii) If f and g have a constant term we assume without loss of generality that p 0 = 0 and we add the following requirement: if p j is an extremal point of conv(p 1 , p 2 , . . . , p t-1 ) then 2p j is not in the support of f and g.

We do not know how to prove a linear upper bound assuming only (i). Condition (ii) is satisfied in particular when the points in Mon(f ) = Mon(g) are convexly independent.

The interesting case, which we consider first, is when f and g have a constant term but f g + 1 has no constant term. As explained above we assume that p 0 corresponds to the constant terms of f and g, i.e., p 0 = 0. Under these hypotheses we have the following result.

Proposition 2. Under assumptions (i) and (ii),

where (p i ) i∈I is the subset of those monomials in Mon(f ) which appear in f g + 1 with a nonzero coefficient.

Proof. We first prove the inclusion from left to right. Since f g + 1 has no constant term, all monomials of f g + 1 are of the form p i + p j where i ≥ 1 or j ≥ 1. Consider first the case where i and j are both nonzero. If i = j this monomial appears in the right-hand side, and if i = j it is the middle point of two points (namely, 2p i and 2p j ) appearing in the right-hand side. The remaining case is when i = 0 or j = 0. If e.g. j = 0 we have p i + p j = p i and we see from the definition of I that this monomial also appears in the right-hand side. Now we prove the inclusion from right to left. Again by definition of I, all the p i with i ∈ I are monomials of f g + 1. Hence it remains to show that conv(2p 1 , . . . , 2p t-1 ) ⊆ Newt(f g + 1).

The left-hand side can be written as conv((2p j ) j∈J ) where the (p j ) j∈J form a convexly independent subset of {p 1 , . . . , p t-1 }. Any monomial of the form 2p j with j ∈ J appears in f g + 1 with a nonzero coefficient because it can be obtained in a unique way by expansion of the product f g. Assume indeed that 2p j = p i + p k with i = k. Then p j is the middle point of p i and p k . If i ≥ 1 and k ≥ 1, this is impossible by construction of J. If i = 0 or k = 0, this is also impossible by hypothesis (ii). We thus have conv((2p j ) j∈J ) ⊆ Newt(f g + 1), and the proof is complete.

We note that this proposition does not hold without assumption (ii), as shown by the following example: take

The monomial X 3 Y 3 is a vertex of Newt(f g + 1) but is not of the form p i or 2p j prescribed by Proposition 2.

Theorem 7. Under the same assumptions (i) and (ii) as above, Newt(f g + 1) has at most t + 1 edges where t denotes the number of monomials of f and g.

Proof. We continue to denote the common support of f and g by {p 0 , . . . , p t-1 }. If 0 does not belong to this support then Newt(f g + 1) is the convex hull of {0} and Newt(f g). Moreover, Newt(f g) = Newt(f ) + Newt(g) = conv(2p 0 , . . . , 2p t-1 ).

If 0 is in the support and f g +1 has a constant term then Newt(f g +1) = Newt(f g) has at most t edges (t and not 2t since f and g have the same support).

In the remaining case (0 is in the support but f g + 1 has no constant term) we need to use hypothesis (ii). This case is treated in Proposition 2. At first sight it seems that Newt(f g + 1) can have up to 2(t -1) vertices, but the list of possible vertices can be shortened by picking a convexly independent subsequence. More precisely, write conv(2p 1 , . . . , 2p t-1 , (p i ) i∈I ) = conv((2p j ) j∈J , (p k ) k∈K ) where J ⊆ {1, . . . , t -1} and K ⊆ I are chosen so that the points in this sequence are convexly independent. By the lemma below, |J ∩ K| ≤ 2. As a result, the number of points in the sequence is

Lemma 3. If p, q, r are 3 distinct nonzero points in the plane then the 6 points p, q, r, 2p, 2q, 2r are not convexly independent. This is clear from a picture and can be proved for instance by considering the 4 points 0, p, q, r. There are two cases.

1. If these 4 points are convexly independent, assume for instance that pq is a diagonal of the quadrangle 0prq. Then the line pq separates 0 from r. As a result, r ∈ conv(p, q, 2r).

2. If the 4 points are not convexly independent, assume for instance that r ∈ conv(0, p, q). In this case, 2r ∈ conv(2p, 2q, r).