published or not. The documents may come

Spatial distribution of mechanical forces and ionic flux

in electro-kinetic instability near a permselective membrane

Spatial distribution of mechanical forces and ionic flux in electro-kinetic instability near a permselective membrane

Pierre Magnico a)

I. INTRODUCTION

Ionic transport in the presence of a charged wall is used in numerous applications such as electro-dialysis, fuel cells, micro-mixers, biochemical analysis, etc. [START_REF]Electrokinetics and Electrodynamics in Microsystems[END_REF][START_REF] Kwak | Enhanced salt removal by unipolar ion conduction in ion concentration polarization desalination[END_REF][START_REF] Nikonenko | Desalination at overlimiting currents: State-of-the-art and perspectives[END_REF][START_REF] Chang | Nanoscale electrokinetics and microvortices: How microhydrodynamics affects nanofluidic ion flux[END_REF] In microand nano-technologies, electro-kinetic mechanisms, such as electro-osmosis and electro-convective instability, play a role all the more important that the specific surface is important and the Debye length to the characteristic length ratio is not so small, i.e., the volume where the electro-neutrality is not fulfilled cannot be neglected.

In membrane systems, ionic separation is characterized by three electro-kinetic regimes depending on the magnitude of the electric potential drop applied to the system. The transport properties such as the ionic conductivity are related to the profile of the non-equilibrium electric double layer (EDL) against the membrane surface. At small potential drop, the electric conductivity is constant (Ohmic regime). As the potential drop increases, the conductivity decreases to a very small value (limiting regime), i.e., the current density reaches an asymptotic value called limiting current. These two regimes have been intensively studied by means of the Nernst-Planck equation coupled to the Poisson equation. [START_REF] Manzanares | Numerical simulation of the nonequilibrium diffuse double layer in ion-exchange membranes[END_REF][START_REF] Zabolotskii | Steady-state ion transport through a three-layered membrane system: A mathematical model allowing for violation of the electroneutrality condition[END_REF] It is observed in particular that in the limiting regime, the EDL splits into two parts: a quasi-equilibrium EDL against the membrane surface and an extended space charge (ESC) layer at the outer edge of a) pierre.magnico@univ-amu.fr the quasi-equilibrium EDL. If the solvent leakage is to be taken into account, the transport of solvent is modeled by means of a linear relation between the solvent velocity and the current density. So the solvent motion as the ion transport is described by a 1D approach. [START_REF] Borges | Investigation of the mass transfer processes during the desalination of water containing phenol and sodium chloride by electrodialysis[END_REF] If the potential drop is high enough, experimental investigations show that the electric conductivity and the current density increase with the potential. This third regime, named over-limiting regime, received several explanations such as water dissociation in the depleted diffuse layer, protonation and deprotonation of the ion-exchange groups at the interface, gravitational effect induced by the ionic concentration gradient or by the joule effect, heterogeneity of the membrane surface (conductivity and rugosity), and electro-kinetic instability along homogeneous membrane. [START_REF] Belashova | Overlimiting mass transfer through cation-exchange membranes modified by Nafion film carbone nanotubes[END_REF][START_REF] Pismenskaya | Evolution with time of hydrophobicity and microrelief of a cation-exchange membrane surface and its impact on overlimiting mass transfer[END_REF][START_REF] Nikonenko | Intensive current transfer in membrane systems: Modelling, mechanisms and application in electrodialysis[END_REF][START_REF] Danielsson | A model for the enhanced water dissociation on monopolar membranes[END_REF][START_REF] Balster | Morphology and microtopology of cationexchange polymers and the origin of the overlimiting current[END_REF][START_REF] Tanaka | Ion exchange membranes fundamentals and applications[END_REF][START_REF] Shapiro | Transport properties of highly ordered heterogeneous ion-exchange membranes[END_REF][START_REF] Belova | Effect of anion-exchange-membrane surface properties on mechanisms of overlimiting mass transfer[END_REF] However, several experimental visualizations have shown that the electro-kinetic instability is the main mechanism at the origin of the over-limiting regime. [START_REF] Kim | Concentration polarization and nonlinear electrokinetic flow near a nanofluidic channel[END_REF][START_REF] Kwak | Shear flow of an electrically charged fluid by ion concentration polarization: Scaling laws for electroconvective vortices[END_REF][START_REF] Rubinstein | Direct observation of a non-equilibrium electro-osmotic instability[END_REF][START_REF] Chang | Understanding electrokinetics at the nanoscale: A perspective[END_REF][START_REF] Kwak | Microscale electrodialysis: Concentration profiling and vortex visualization[END_REF] Electro-convective instability was first explained by the presence of the ESC region and a heterogeneous electric field (spherical ion-exchange resin and heterogeneous ionexchange membrane). [START_REF] Mishchuk | Electro-osmosis of the second kind near the heterogeneous ion-exchange membrane[END_REF][START_REF] Rubinstein | Electroconvection at an electrically inhomogeneous permselective interface[END_REF][START_REF] Dukhin | Electrokinetic phenomena of the second kind and their applications[END_REF] However, in the case of homogeneous membrane, two theoretical approaches carried out by Rubinstein and Zaltzman [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF][START_REF] Rubinstein | Electroconvective instability in concentration polarization and non-equilibrium electro-osmotic slip[END_REF][START_REF] Rubinstein | Electro-osmotic slip of the second kind and instability in concentration polarization at electrodialysis membranes[END_REF][START_REF] Zaltzman | Electro-osmotic slip and electro-convective instability[END_REF] have proved the existence of the electro-kinetic instability. An asymptotic analysis for small Debye lengths leads to the determination of the expression of the velocity slip at the outer edge of the ESC region. A linear instability analysis was also developed to estimate the potential drop threshold and the vortex size. Suh 28 This paper is devoted to the numerical investigation of electro-kinetic instability in a polarization layer next to a cation-exchange membrane. An analysis of some properties of the electro-kinetic instability is followed by a detailed description of the fluid flow structure and of the spatial distribution of the ionic flux. In this aim, the Stokes-Poisson-Nernst-Planck equation set is solved until the Debye length scale. The results show that the potential threshold of the marginal instability and the current density depend on the logarithm of the concentration at the membrane surface. The size of the stable vortices seems to be an increasing function of the potential drop. The fluid motion is induced by the electric force along the maximum concentration in the extended space charge (ESC) region and by the pressure force in the region around the inner edge of the ESC layer. Two spots of kinetic energy are located in the ESC region and between the vortices. The cationic motion, controlled by the electric field and the convection, presents counter-rotating vortices in the stagnation zone located in the fluid ejection region. The anion transport is also characterized by two independent layers which contain counter-rotating vortices. The first one is in contact with the stationary reservoir. In the second layer against the membrane, the convection, and the electric field control, the transversal motion, the Fickian diffusion, and the convection are dominant in the longitudinal direction. Finally, the longitudinal disequilibrium of potential and pressure along the membrane is analyzed. performed a thorough instability analysis over a wide range of several parameters from the monotonic to the oscillatory mode.

Even if these theoretical studies shed a new light on electro-convective instability in membrane systems, it is necessary to solve the fluid and ionic transport equations in order to reach all the details of the electro-convective instability and to expand the domain of validity of the investigation. In this aim, Demekhin research group [START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF][START_REF] Chang | Competition between Dukhin's and Rubinstein's electrokinetic modes[END_REF][START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaos motion[END_REF][START_REF] Demekhin | Thermoelectroconvection near charge-selective surfaces[END_REF][START_REF] Demekhin | Three-dimensional coherent structures of electrokinetic instability[END_REF] studied the temporal evolution of the instability expansion (self-similarity and ESC layer deformation), the instability mode interaction in the presence of rugosity, the instability regimes (from the bifurcation analysis for steady instability to the chaotic motion), the 3D electro-convective pattern in the function of the potential drop. In the same spirit, Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF][START_REF] Pham | Helical vortex in three-dimensional electrochemical systems with ion-selective membranes[END_REF] analyzed the current hysteresis at the transition limiting/over-limiting regime and helical hydrodynamic motion. Druzgalski et al. [START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF][START_REF] Druzgalski | Statistical analysis of electroconvection near an ion-selective membrane in the highly chaotic regime[END_REF] studied the chaotic mixing. The authors used ensemble averaged equations to analyze the spatial distribution of the energy dissipation, the energy spectrum, and the fluctuating forces. Direct numerical simulations are also used in the chemical engineering context. Kwak et al. [START_REF] Kwak | Enhanced salt removal by unipolar ion conduction in ion concentration polarization desalination[END_REF][START_REF] Kwak | Shear flow of an electrically charged fluid by ion concentration polarization: Scaling laws for electroconvective vortices[END_REF] and Urtenov et al. [START_REF] Urtenov | Basic mathematic model of overlimiting transfer enhanced by electroconvection in flow-through electrodialysis membrane cells[END_REF] studied the electro-convective instability in the case of a pressure driven shear flow between two ion-exchange membranes. A scaling of the vortex size and of the vortex advection speed was found, and several modes of instability have been investigated by chrono-potentiometry from the Dukhin-Mishchuk mode to the unsteady Rubinstein-Zaltzman one.

So, after a long period of time spent on the linear stability analysis, direct numerical simulations have given an insight on the hydrodynamic structure and on the ionic flux contributions. However, these results describe the instability evolution in the transient and chaotic mode and also under pressure driven flow. To my knowledge, results dealing with marginal instability computed by direct numerical simulations are published by Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] and Demekhin et al. [START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaos motion[END_REF][START_REF] Demekhin | Three-dimensional coherent structures of electrokinetic instability[END_REF] only. So, the goal of the present work is to describe, in more detail and in the simplest bi-dimensional situation, the contribution of each term of the transport of momentum and ionic species at steady-state.

After a validation step, the electro-convective instability is studied at the global level. The coherence between the noflux of co-ions and the counter-ion imposed concentration at the membrane surface is addressed first. The sensitivity of the current density and of the potential threshold toward the counter-ion imposed concentration is shown and a correlation is proposed. The simulations also show that the range of stable roll size depends on the potential drop value. In a second step, a detailed analysis of the contribution of each term of the transport equations is performed. The spatial distribution of the electric and of the pressure force is compared first to the spatial distribution of the electric charge and of the kinetic energy. We show how they act on the fluid motion and at which location they are dominant. Then the spatial analysis of each contribution of the ionic flux (diffusion, convection, and electro-migration) is carried out to explain the motion of each ionic species. The results show in particular the presence of counter-rotating vortices of species flux even for the counterion despite the high electric force directed to the membrane surface. Finally, the longitudinal disequilibrium of pressure and potential, leading to the electro-kinetic instability, is studied at the roll level. We show at the same time the agreement between the macroscopic analysis and the microscopic one developed previously.

II. THEORY

A. Model and governing equations

In this work, we use the two-dimensional model system described in the studies of Druzgalski et al. [START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] and Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] (Fig. 1). An electrolytic solution lies between a cationexchange membrane and a bulk space at rest. In the bulk, located at L x from the membrane, the ionic concentration c o is homogeneous. The longitudinal length of the domain is L z . At the interface, the cationic concentration is constant and equal to c + interf and a no-flux of anion is imposed. A potential drop ∆Φ is applied through the domain in the transversal direction x. Periodic condition is imposed in the longitudinal direction z.

The incompressible fluid motion is described by the Stokes equation [Eq. ( 1)] coupled with the continuity equation [Eq. (2)]. The ionic transport is governed by the Nernst-Planck equation [Eq. ( 3)] and the electric potential is related to the local ionic concentration by the Poisson equation [Eq. ( 4)],

0 = -∇P + µ∆ U -ρ f z + c + + z -c -∇Φ, (1)
∇ • U = 0, (2)
d dt c ± = ∂ ∂t c ± + U • ∇c ± = -D ± ∇ ∇c ± + z ± c ± k B T ∇Φ , (3)
-ε o ε r ∆Φ = Q = z + c + + z -c -, ( 4
)
where P is the pressure, µ is the dynamic viscosity, U is the fluid velocity, ρ f is the fluid density, z ± is the charge of the cation and of the anion, c ± is the concentration of the cation and of the anion, Φ is the electric potential, t is the time, D ± is the diffusion coefficient of the cation and of the anion, Q is the electric charge density, ε o and ε r are the permittivity of the vacuum and the relative permittivity of the solution, respectively, and k B and T are the Boltzmann constant and the temperature, respectively. In Eq. ( 3), the interaction between the current density and the induced magnetic field is neglected. The speed of light is far much higher than the diffusional velocity computed at the domain scale. Therefore, we assume that the magnetic potential is constant in time. These equations are scaled as in the work of Druzgalski et al. [START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] The characteristic parameters are the transversal length L x , the diffusional time t diff = L 2

x D, the diffusional velocity U diff = D/L x , the thermodynamic potential Φ T = k B T e, the "osmotic" pressure P o = µD L 2

x , and the bulk concentration c o . Here e is the elementary charge. In the following, we assume that D + = D = D so that for a symmetric univalent electrolyte, Eqs. ( 1)-( 4) become in their dimensionless form

0 = -∇P + ∆ U - Pe 2ν 2 c + -c -∇Φ, (5)
∇ • U = 0, ( 6
)
d dt c ± = ∂ ∂t c ± + U • ∇c ± = -D ± ∇ ∇c ± ± c ± k B T ∇Φ , ( 7
)
-2ν 2 ∆Φ = Q = c + -c -, (8)
where Pe and ν are the Peclet number and the dimensionless Debye length, respectively,

Pe = ε o ε r µD k B T e 2 , (9)
ν = λ d L x with λ d = ε o ε r k B T 2e 2 c o 1/2 . ( 10
)
For simple ions, the value of Pe and D is around 0.5 and 10 9 m 2 /s, respectively. In this work, besides these two parameter values, ν is imposed equal to 10 3 . So, if L x = 0.1 mm, then λ d = 100 nm, c o = 10 5 M (T = 300 K), and the diffusional time is 10 s. From now, we will use the dimensionless variables. Therefore, the value of dimensionless transversal length L x and of the dimensionless bulk concentration c o is equal to 1.

The boundary conditions used at the fluid/membrane interface are the no-slip condition, an imposed value of the cationic concentration c + interf , a no-flux for the anion, and a null electric potential. At the bulk boundary, no-slip condition is imposed and the potential value ∆Φ is applied. In the longitudinal direction, the periodic condition is used. In this study, c + interf , L z , and ∆Φ vary from 1 to 32, from 6/7 to 6, and from 0 to 50 (0-1.3 V), respectively. In the following, the ionic flux and the current density are defined in their dimensionless form,

J ± = c ± U -∇c ± ± c ± ∇Φ
and

I = I x Ω = 1 A Ω J + x -J - x dA, ( 11
)
where A is the area of the computational domain Ω. The expression of I is the definition used in the work of Druzgalski et al. 36

B. Numerical method

As in the work of Pham et al., [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] the collocation finite volume is used. The Stokes and the continuity equations are solved by means of Darwish's method. [START_REF] Darwish | A coupled finite volume solver for the solution of incompressible flows on unstructured grids[END_REF] This method consists in the use of the Rhie-Chow interpolation in the continuity equation. This allows us to solve the Stokes and the continuity equations at the collocation points. Darwish's method also increases the robustness of the resolution because contrary to other algorithms like SIMPLE, the velocity at the surface of the control volume is computed implicitly in the divergence free equation without estimating a pressure correction. At the membrane surface and at the bulk boundary, the pressure is linearly extrapolated. The Poisson-Nernst-Planck (PNP) equation set is solved by means of the Newton-Raphson method in fully implicit form. The two equation sets are solved alternatively. The second order space discretization is used. In the Nernst-Planck equation, the convective term is computed with the Quick scheme. As regards the time discretization, the backward second order scheme is used. However, unlike Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] a structured grid is employed instead of a triangular mesh.

Owing to the high gradient near the interface, the mesh must be refined in the x direction. In this aim, a geometric series is used. This refinement induces an increase of the aspect ratio (longitudinal stretching) of the control volume as we move closer to the membrane surface. In order to avoid this problem and to use cells with an aspect ratio value between 0.75 and 2, the control volumes are divided up into two equal parts in the direction z when the aspect ratio is greater than a desired threshold value (Fig. 2). The number of longitudinal refinement levels varies from 4 to 6 according to the value of ν and to the needed accuracy. At the boundary between two refinement levels, a second order interpolation is used.

In the case of simulations carried out with c + interf ≤ 4, ν = 10 3 , and L z = 6, 180 and 384 nodes are used in the transversal and longitudinal directions, respectively. In the transversal direction, the common ratio of the geometric series is 0.0234. In the longitudinal direction, the refinement threshold is 2. FIG. 2. Sketch of the numerical mesh. Refinement level from N to N + 2.

In the case of simulations carried out with 8 ≤ c + interf ≤ 32, ν = 10 3 , and L z = 6, 250 and 485 nodes are used in the transversal and longitudinal directions, respectively. The common ratio of the geometric series is 0.02. In the longitudinal direction, the refinement threshold also is 2. For these two series of simulation, Eqs. ( 5)-( 8) are integrated in time with a time step of 10 4 .

The mesh independency has been verified for ∆Φ ≤ 40 as concerns the current density. The number of nodes in the x-direction and the value of the geometric series are chosen in order to have at least 10 nodes in the quasi-equilibrium EDL with a small variation rate of the cell size. In the z-direction, the node number is chosen in order to have a cell aspect ratio equal to 0.75 at the bulk boundary. Numerical simulations with a refinement threshold of 2 and 2.5 give the same results. As regards the time step, the value is not small enough for accurate simulations. In this work, the transient phenomena are not important. In the unsteady instability region, the simulations are performed in order to have mean values. As concerns the self-similarity validation, the time step values of 10 5 and 10 4 give the same results. However, in order to verify the hysteretic behavior with c + interf = 1, the common ratio of the geometric series and the time step are decreased to 0.005 and 10 5 , respectively.

C. General remarks about the electro-membrane system conductivity vs ∆Φ

Before discussing the numerical results, the main transport properties of membrane systems must be recalled first. Figure 1 displays the charge density and the streamlines induced by a steady-state electro-kinetic instability when ∆Φ = 20 (c + interf = 2, ν = 10 3 , and L z = 6). The electric charge density is positive in the whole domain but for the sake of clarity, the color map lies between 10 3 and 3 × 10 3 . Two regions characterize the spatial distribution of the ionic species: the diffusion layer, where the electro-neutrality is fulfilled (x ∈ [0.15, 1]), and the non-equilibrium electric double layer (EDL) of width 0.15 against the fluid/membrane interface. This region of high electric charge density has a much larger width than the Debye length, i.e., the width of the EDL at a charged interface in equilibrium with an ionic solution. The non-equilibrium EDL is usually divided into 2 sub-layers: the extended space charge (ESC) layer against the diffusion layer and the quasiequilibrium EDL between the ESC layer and the membrane surface. The form of the polarization layer, i.e., the presence and the width of the non-equilibrium EDL, and the conductivity properties also depend on the potential drop value (see Fig. 2 in Ref. 24, for example).

In the Ohmic regime (region I in Fig. 3), the ESC region does not exist and the diffusion layer, where the concentration gradient is constant, takes place in most of the polarization layers. The concentration gradient increases with ∆Φ and the electric conductivity (dI/d (∆Φ)), depending on the ionic concentration, decreases with ∆Φ. At the same time, the quasi-equilibrium EDL width increases. When ∆Φ reaches a first threshold value, the ionic concentration vanishes at the fluid/membrane interface and the limiting regime is reached (region II in Fig. 3). Assuming a null value of c + and a no-flux FIG. 3. Current density versus electric potential drop. Comparison between numerical results obtained in this work and those published by Druzgalski et al. [START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] and Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] Black curve: computations carried out in the case of a quiescent solution. Four regimes of electric conductivity: (I) Ohmic regime, (II) limiting current regime, (III) steady-state over-limiting regime, (IV) unsteady over-limiting regime.

of anions at the membrane interface, and also assuming the electro-neutrality, we obtain from Eq. ( 11) a simple expression of the limiting current density,

I lim = 2D + c o L x . (12)
The dimensionless expression of I lim has a value of 2 as shown in Fig. 3.

As ∆Φ increases, the concentration gradient in the diffusion layer continues to increase. The point, where the ionic concentration vanishes, moves away from the membrane surface. Between this point and the outer edge of the quasi-equilibrium EDL takes place the ESC layer. This layer comes from a part of the quasi-equilibrium EDL whose width decreases. The presence of a maximum of electric charge density and of a minimum located at the outer edge of the quasi-equilibrium EDL characterizes the ESC layer. As long as ∆Φ is smaller than a second threshold value, the width of the non-equilibrium EDL increases and the ionic solution remains at rest. This means that the ionic transport can be described by unidimensional considerations. As we go through this second threshold (∆Φ * ), the apparent conductivity increases suddenly with ∆Φ (region III in Fig. 3). Beyond this threshold, the electro-convective instability is induced by the presence of longitudinal heterogeneity in the non-equilibrium EDL region. The ionic and the fluid transport become multidimensional but steady. By increasing ∆Φ, a third threshold is reached and the ion transfer becomes unsteady (region IV in Fig. 3).

We must notice that close to the thermodynamic equilibrium (∆Φ = 0), I is negative because the Donnan potential is not taken into account. This approximation is discussed in Sec. III B.

III. RESULTS

A. Validation

In order to validate the numerical method, results are compared to those published by Druzgalski et al., [START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] Pham et al., [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] and Demekhin et al. [START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF] In Refs. 34 and 36, the two-dimensional membrane system is identical to the system described here with the same boundary condition. However, in the work of Pham et al., [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] a small rugosity takes place at the fluid/membrane interface. The published numerical results are obtained with c + interf = 2, ν = 10 3 , and Pe = 0.5. In the work of Demekhin et al., [START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF] the polarization layer is located between two cationexchange membranes with c + interf = 5, ν = 5 × 10 -4 / √ 2, and Pe = 0.25. Figure 3 shows the agreement between the published data [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF][START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] and those computed in this work. In this figure, the black solid curve in the Ohmic regime and the limiting one is computed with a 1D simulation.

Let us consider the results of Druzgalski et al. [START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] The domain length value used by the authors is 6. Here, the currentpotential curve is computed first at ∆Φ = 20. When the stationary solution is reached, i.e., J -

x Ω ∼ 0, the potential drop is increased to a higher value until ∆Φ = 42.5. The value of L z is 6. For ∆Φ = 20, a first 1D solution is computed with the Stokes-PNP equation set and then reported over all the computational domains. So at the first time step, all the variables are independent of z. In order to accelerate the instability development, a small amplitude white noise is superimposed on the concentration and the electric potential over the computational domain (c

± (z, x) = c ± (x) × (1 + a × Rand) and Φ(z, x) = Φ(x) × (1 + a × Rand),
where Rand is a random number whose value lies between 0 and 1 and a is the perturbation amplitude). As in Ref. 36, we obtain a vortex pair size of 6/5. In the steady instability regime, the roll size does not change with ∆Φ. If ∆Φ > 42.5, the vortices are no more stable and the mean transversal current density displays a non-periodic temporal variation. So, in this unsteady regime, the current density is averaged over the time and the computational domain. The vertical bars represent the root mean square of the signal. A second set of computation is carried out from ∆Φ = 25 in the decreasing direction with L z = 1.2 = 6/5 in order to find the potential threshold of marginal instability corresponding to the vortex pair size.

A 2nd test consists in the comparison with the results of Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] The authors use a domain length value of 1. In the present work, the profiles obtained with L z = 1.2 are used as the initial condition. A homothetic transformation is applied in the z direction in order to obtain L z = 1 with the same number of computation nodes. The authors [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] find a hysteresis of current density at the transition to the overlimiting regime. In the increasing potential drop direction, the threshold value lies between 20 and 21 and the threshold one of 19.1 is found in the decreasing direction. In the present work, a subcritical bifurcation is also observed but the hysteresis amplitude is much smaller. The profiles computed at ∆Φ = 19 are used as the initial condition. Switching ∆Φ to a value of 19.25, U 35. This small difference in potential drop threshold of 0.2 is in accordance with the value found by Demekhin et al. [START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaos motion[END_REF] A third test deals with the self-similar behavior of the current-voltage curve. The test is carried out with c + interf = 1, ν = 10 3 , and ∆Φ = 20, 30, 40, 50, i.e., from the marginal instability threshold to the unsteady instability one. Demekhin et al. [START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF] show that at intermediate time before the instability development, the diffusion layer and the current density follow a self-similar evolution. The authors rescale the variables with the diffusional length √ 4Dt. The new dimensionless Debye length ν , current density I , and electric potential drop ∆Φ are defined as

ν = ν/2 √ t, I (t) = 2I (t) √ t, ∆Φ (t) = ∆Φ -I (t)/2. ( 13
)
To verify the self-similar evolution through the currentvoltage curve I (t) vs ν ∆Φ , the procedure described in Ref. 29 is not followed exactly. In this article, [START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF] the initial condition consists in 1D profiles computed at ∆Φ = 0 to which is superimposed a white-noise perturbation in concentration. Then the electric potential drop is switched to a prescribed value lying between 50 and 200. In the present work, the boundary condition c + interf = 1 is used because the condition ∆Φ = 0 is not consistent with the thermodynamic equilibrium between the membrane and the electrolytic solution (i.e., no ionic flux) if c + interf 1. With c + interf = 1, the initial condition is a flat profile of concentration and potential, i.e., the equilibrium EDL is not present against the membrane (see comments in Sec. III B).

Figure 4 shows the variation of I (t) vs ν I (t)/2 instead of I (t) vs ν ∆Φ . [START_REF] Demekhin | Linear and nonlinear evolution and diffusion layer selection in electrokinetic instability[END_REF] The circles stand for the universal current-voltage characteristics and the solid lines represent the numerical results for the four values of the potential drop. The characteristics have the following expression:

I (t) = 2.16 + 3.14 ν I (t)/2, ( 14
)
where the value 2. increases exponentially (region III). Finally, the fourth one (t > 0.2 or ν I (t)/2 < 0.016), the mean kinetic energy value is is greater than 1, and the ionic transport is no more one-dimensional (region IV).

B. Electro-convective instability property effect of c + interf and ∆Φ

In this section, two aspects are discussed: the dependence of the instability threshold on the cation concentration at the membrane surface and the dependence of the roll size on the imposed potential drop.

As regards the first point, until now the value of c + interf was 1 or 2. However, the selectivity of a membrane system is ideal if the membrane charge, i.e., c + interf , is much greater than the ion concentration in the bulk (c o ). In the literature, the imposed value of c + interf is equal to 1 or 2 in Refs. 34-38, and the value is equal to 4 or 5 in Refs. 29-33. Rubinstein research group found, by means of an asymptotic expansion, that the electro-kinetic of the second kind does not depend on this parameter. But Fig. 5 shows that the potential threshold (∆Φ * ) and the current density are dependent on c + interf . The present investigation is carried out with c + interf = 1, 2, 4, 8, 16, 32 for two domain lengths L z = 6 5 and L z = 6/7. The length is chosen to impose the vortex pair size. In the Ohmic and the limiting regime, the curves are computed with c + interf = 1, 2, 8. In the Ohmic regime, the figure displays a negative value of the current density at ∆Φ = 0. This comes from the inconsistency between the condition ∆Φ = 0 and the null current density due to the Donnan potential induced by the presence of the fixed charge density inside the membrane. This potential is computed by assuming (1) the thermodynamic equilibrium between the membrane and the ionic fluid and (2) the electroneutrality inside the charged membrane. The thermodynamic equilibrium is defined by the equality of the chemical potential over the domain µ ± (x = 0) = µ ± (x) = Ln(c ± ) + z ± Φ. From these assumptions and definition, the Donnan potential has the following expression:

Φ d = Ln -Y + √ Y 2 + 4 2 , ( 15
)
where Y is the dimensionless concentration of fixed charge in the membrane.

Φ
(Φ interf = -Ln c + interf ).
Assuming that the anionic flux J is null, the current density is equal to the cationic flux. J = 0 means that the anions are in thermodynamic equilibrium. Using the equality of the anionic chemical potential at the interface (x = 1) and at the reservoir edge (x = 0) and using the boundary condition on the potential, one obtains c - interf = c o = 1. Let us consider the cations. The current density is defined by means of the chemical potential gradient

I = J + = -c + ∂/∂x µ + .
Integrating the current over the polarization layer, we find

I × 1 ∫ 0 1/c + (x)dx = -µ + (x = 1) = -Ln c + interf . Therefore, the current density is negative if c + interf > 1 and vanishes if c + interf = 1.
To conclude, the thermodynamic equilibrium is fulfilled (J ± = 0) if the potentials at infinity on both sides of the membrane are equal. [START_REF] Manzanares | Numerical simulation of the nonequilibrium diffuse double layer in ion-exchange membranes[END_REF] Or we must use the boundary condition c + interf = 1 if one reservoir is considered. However this problem has no influence if the imposed potential drop is higher than Φ d , i.e., the operating point is in the limiting and over-limiting current regime as shown when ∆Φ > 10.

Figure 5 shows that the potential threshold ∆Φ * is an increasing function of c + interf and a decreasing function of L z , i.e., the vortex pair size so that at a fixed value of ∆Φ, the current density is a decreasing function of c + interf and an increasing function of L z . We can also observe that as ∆Φ increases, the current density is less sensitive to c + interf so that close to the unsteady instability threshold, I is independent of c + interf . Figure 6(a) shows that an empirical relation in the form ∆Φ * = a + b × Ln c + interf agrees the numerical results. Figure 6(b) shows that at a fixed value of ∆Φ, the current density follows the same logarithmic variation with c + interf . The values of a and b are gathered in Table I. From Fig. 5, we can make another remark: the width of the potential window, where the steady-state instability is a solution of the problem, decreases with the roll size and the marginal instability does not exist when ∆Φ > 35 if the vortex pair size is 6/7.

The logarithmic behavior of the potential threshold can be compared to its dependence on the length ratio ν. [START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaos motion[END_REF] In Fig. 5, the authors [START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaos motion[END_REF] show by means of a linear stability analysis that in the range of ν lying between 10 4 and 10 3 , the potential threshold decreases logarithmically from 22 to 19, confirming the analysis of Rubinstein and Zaltzman. [START_REF] Zaltzman | Electro-osmotic slip and electro-convective instability[END_REF] The range of variation of ∆Φ * is exactly the same as the range shown in the present work. So even if c + interf varies over a range smaller than two orders of magnitude, the dependence of ∆Φ * on c + interf must be taken into account.

In the following, the influence of ∆Φ on steady periodic solutions is addressed. A thorough study of stability was performed by Demekhin's research group [START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaos motion[END_REF] in the case of a polarization layer between two cationic exchange membranes. In particular, the authors performed a stability analysis of steady space periodic solutions by means of a small amplitude white-noise perturbation in order to investigate the transition from the steady solution to the time periodic one and the chaotic mixing as a function of ∆Φ in the case of supercriticality (Pe = 0.1).

TABLE I. Parameter values of the equation

F = a + b × Ln c + interf .
Electric potential threshold (F = ∆Φ * ) for two domain lengths L z . In the present work, the influence of ∆Φ on steady periodic solutions is investigated in the case of the subcritical bifurcation (Pe = 0.5). Unlike the work of Demekhin et al., [START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaos motion[END_REF] unsteady solutions are not considered. In fact, if ∆Φ = 20 with c + interf = 2, ν = 10 3 , and L z = 6, the vortex pair size of the steady-state solution is L z /5, while the size is L z /6 if ∆Φ = 30 with c + interf = 8, ν = 10 3 , and L z = 6. So it seems that the steady roll size and the size range are increasing functions of ∆Φ. In order to verify these observations, an investigation is carried out close to ∆Φ * (∆Φ = 20) and far from ∆Φ * (∆Φ = 30) with c + interf = 2, ν = 10 3 , and L z = 6. Figure 7 displays the temporal variation of the amplitude A(λ) of six wavelengths (λ = L z /3, L z 5, L z /6, L z /7, L z /8) during the development of the instability when ∆Φ = 20 and 30. The main wavelengths are chosen in order to simplify the investigation of the transitory process in the linear domain. These curves are obtained from the FFT analysis of U z . Initially a signal composed of these six sinusoidal waves in the z direction is added to the 1D solution at ∆Φ = 20 and 30. The amplitude of the signal is 0.1%. The perturbation concerns c ± and Φ. The simulations are carried out in the linear domain, i.e., during the exponential increase of U

λ o L z /4 L z /5 L z /5 L z /6 λ o L z /5 L z /4 L z /6 L z /5 λ final L z /4 L z /5 L z /5 L z /5
wave length decreases. So the wavelengths, the amplitude of which increase the most rapidly, decrease as ∆Φ increases.

In

C. Detailed description of the hydrodynamic

The spatial distribution of the mechanical forces and of the local coupling between the forces and the velocity field has not been described partially in the past. [START_REF] Zaltzman | Electro-osmotic slip and electro-convective instability[END_REF][START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] In Ref. 27, the authors used an asymptotic expansion to explain the instability mechanism. In the discussion, the transversal profile of the longitudinal component of the two forces shows that the pressure is the main contribution of the instability, far ahead the electric force. On contrary, in Ref. 34, the two forces act on the instability mechanism. In the discussion, a 1D approach is also used as regards the longitudinal component. In the following, we use the results computed with the parameter values ∆Φ = 25, c + interf = 1, L z = 6/5, and ν = 10 3 . Computations carried out with other values show similar properties. Half of the domain (z ≤ 0.6) is taken into account in order to simplify the comments owing to the symmetry of the problem.

Figure 8 shows the spatial distribution of the charge density (Q) and of the kinetic energy (E c ) on which the streamlines are superimposed. The vortex centers are located at (x = 0. increases quickly in the region z < 0.3 and reaches a maximum at z = 0.45. As the ESC layer expands, its mean charge and E c decrease. So the region of high fluid acceleration is correlated to the presence of a maximum charge in the ESC layer. The presence of a secondary maximum of E c is explained by the flow convergence to the symmetry axis when x < 0.11. As soon as the streamlines diverge, i.e., when x > 0.11, U

x decreases as x increases. For the same reason, we should see a spot of kinetic energy in the region around z = 0. However, the distance between two neighboring vortices in this region is higher than around z = 0.6. So the convergence of the streamlines has less effect.

Figure 9 displays the difference in magnitude between the electric force and the pressure force in the longitudinal

direction ( F E z -F Pr z = -c + -c -∇Φ z --∇P z
) and in the transversal one. The Stokes equation describes the equilibrium between the viscous, the electric, and the pressure force. So, the figure shows where one of the two forces is dominant and is in equilibrium with the viscous force. The figure shows two layers along the fluid/membrane interface and two minima on both sides of the symmetry axis.

The first layer with a positive value (preponderance of the longitudinal electric force) is at the same transversal position than the high kinetic energy layer. Along the membrane from z = 0.1 to z = 0.5, this layer follows the displacement of the maximum of charge density Q in the ESC region, i.e., the distance from the membrane surface increases with z. F E z is positive of course. The pressure increases with z, reaches a maximum value at z = 0.5, and then decreases around the symmetry axis. F Pr z is therefore negative and its magnitude increases with z until z = 0.5: F Pr z acts against the fluid motion. In the second layer, F Pr z is dominant and

F Pr z > 0.
However F E z is also positive all along the membrane (z < 0.6). Its magnitude is smaller than F Pr z owing to the boundary condition imposed on Φ at the membrane surface. F Pr z and F E z increase with z, reach a maximum at z ∼ 0.45, and decrease to 0 at z ∼ 0.6. Along the membrane, the width of the 2nd layer increases with the nonequilibrium one and in particular with the transversal position of the minimum of Q.

At a distance greater than 0.075 from the interface, Q is very small and F E z is no more dominant. So the fluid motion, induced by the electric force at the bottom of the vortices and by the friction in the transversal direction, produces a negative longitudinal pressure force. At x = 0.1, as long as the amplitude of F Pr z remains small (z < 0.35), U z increases. Then F Pr z decreases quickly inducing a high decrease of U z as z increases. F Pr z reaches a minimum in the region z ∼ 0.5 owing to the symmetry condition at z = 0.6. Around the symmetry axis, the pressure reaches a maximum value if x > 0.05 and the extremum forms a plateau, i.e., it does not depend on z, if x ∈ [0.2, 0.45]. The width of the plateau increases with x from 0.05 to 0.1. At the same time, U x follows a quadratic function of z over a longitudinal range larger than the plateau width. Now let us consider the transversal contribution of the two forces. In the whole domain, F E x is negative which is not the case for F Pr x . Here also, if x < 0.15, we can observe three distinct layers parallel to the fluid/membrane interface: a first layer against the membrane surface where the electric force is dominant, a second one (x < 0.035) where the pressure force is dominant, and a third one (x < 0.15) where the electric force is dominant again. The first layer corresponds to the quasi-equilibrium EDL, the second one is located in the region around the minimum charge density in the ESC layer, and the last one lies in the remaining ESC region. The two outer edges of the 2nd and the 3rd layers, i.e., x = 0.035 and x = 0.15, also correspond to the width of the layers where F Pr z > 0 and and F E x < 0. The difference F E x -F Pr x seems to be large but it is small compared to each force magnitude which is equal to 10 6 in the ESC layer and increases to 10 8 in the quasi-equilibrium EDL owing to the high increase of Q. So we may assume the equilibrium between the two forces in the non-equilibrium EDL, i.e., where U x is close to zero (see the brief discussion in Sec. III E).

F E z > 0,
At the domain center z ∼ 0.6, we observe that the width of the 1st layer increases and the 2nd layer is moved away from the membrane surface and its width also increases. But the 3rd one disappears. The expansion of the first two layers comes from the enlargement of the quasi-equilibrium EDL and the width of the 2nd layer is also increased by the streamlines convergence. The decrease of Q in the ESC region induces a decrease of F E x and explains the absence of the 3rd layer.

D. Comparison with asymptotic result

In Sec. III B, the numerical results show that the potential threshold and the current density depend on c + interf . However, the asymptotic expansion [START_REF] Rubinstein | Electro-osmotic slip of the second kind and instability in concentration polarization at electrodialysis membranes[END_REF][START_REF] Zaltzman | Electro-osmotic slip and electro-convective instability[END_REF] performed in the over-limiting regime leads to a description of the non-equilibrium EDL independent of c + interf over all the current regimes but in the over-limiting one. So, in order to identify the disagreement, the asymptotic results are compared to the numerical ones.

In the case of standard ion concentration and transversal length, the non-equilibrium EDL width is very small (ν = 10 5 if c o = 10 1 M and L x = 0.1 mm, for example). In order to simplify the problem of ion transport and to avoid taking into account the non-equilibrium EDL, Rubinstein and Zaltzman [START_REF] Rubinstein | Electro-osmotic slip of the second kind and instability in concentration polarization at electrodialysis membranes[END_REF][START_REF] Zaltzman | Electro-osmotic slip and electro-convective instability[END_REF] consider the transport of a quasi-electroneutral solution in the diffusion layer with a slip boundary condition at the membrane surface. By matching the asymptotic expansion in the diffusion layer (named "outer" electroneutral layer) and in the non-equilibrium EDL (named "inner" layer), the authors find an expression of the position of the limit between the two layers (x o ) and of the slip velocity (U s ). Therefore, in the simplified model, the fluid/membrane interface is displaced to the outer edge of the ESC layer.

At this step, we must notice that in the Rubinstein notation, the Debye length λ d and the characteristic velocity U o have not the same expression as in the present work,

ν = k b T ε o ε r e 2 c o = √ 2ν
and

U o = ε o ε r µD k b T e 2 U diff = 0.464U diff . ( 16
)
The expressions of x o and of U s are determined by means of the Stokes equation and the asymptotic development of the electric field,

x o = ν 2/3 ∆Φ + 2 3 Ln ν 3 2 √ 2I 2/3 , ( 17
)
U s | x=x o = - 1 9ν 2 dI dz x 3 o , ( 18
)
TABLE III. For three longitudinal positions z are gathered the positions of the minimum of the electric charge density estimated from Eq. ( 19) (x o ) and of the minimum computed in the present work (x c + min ), the position of the outer edge of the ESC layer estimated from Eq. ( 17) (x o ), and the ratio of the numerical velocity to the slip velocity estimated from Eq. ( 18 where I = I(z) is the local current density. In order to test the agreement between these two analytical expressions and the numerical results, the value of x o is determined by computing I = I(x, z) with the numerical data and by matching the value of x with the expression of x o and then the theoretical slip velocity can be compared to the numerical one.

Table III

x o = 2ν 1 c + min - 1 c + interf , ( 19
)
where c + min is the minimum value of the cationic concentration at the inner edge of the ESC layer. In most of the non-equilibrium EDL all along the membrane, the numerical results confirm Rubinstein's hypothesis as regards the mechanical equilibrium in the transversal direction. Neglecting the viscous force in the Stokes equation and using the Poisson equation in order to substitute Q with the Laplacian of the potential, the integration of the Stokes equation leads to

P (x, z) = Pe 2 ∂Φ (x, z) ∂x 2 + b. ( 20
)
This equality is fulfilled if x < x with x and b increasing from 0.05 to 0.12 and from 1730 to 1140, respectively, as z increases from 0 to 0.6. The thermodynamic equilibrium hypothesis is approximately fulfilled in the quasi-equilibrium EDL when c + (x) > 0.05. the region x < 0.023 which corresponds to the layer where

F Pr z >> F E z
. It also corresponds to the region where the vorticity magnitude is the highest. Assuming a Poiseuille profile, 2A × B should be equal to the sum of the electric and the pressure force. Over the three positions, the mean value of 2A × B is around 1.1 × 10 5 , whereas the mean force is equal to 8.5 × 10 4 and 1.1 × 10 5 at x = 0.02 and at the interface, respectively.

The discrepancy lies in the too high value of the estimated slip velocity compared to the numerical ones (see Table III The simplified transport equation of electroneutral species is coupled with the two modified boundary conditions at the membrane surface: the slip condition and the null concentration. Therefore the expression of U s must depend explicitly on c + interf . In the over-limiting regime, the authors [START_REF] Zaltzman | Electro-osmotic slip and electro-convective instability[END_REF] simplified the integral expression of U s by keeping the leading terms. The c + interf dependency may be achieved by a higher order development. However, the dependence of U s on dI/dz does not seem right enough.

E. Description of the cationic transport

The ionic flux depends on three contributions: the Fickian diffusion J ± Fi = -∇c ± , the electro-migration J ± Em = ∓c ± ∇Φ, and the convection J ± Co = c ± U. The transversal cation flux has been briefly discussed by Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] The authors compare three flux profiles along the longitudinal direction close to the membrane surface (x = 0.02) in the presence and in the absence of vortices. They show that the electro-migration is dominant in the two cases. In the presence of vortices, the transversal flux reaches a maximum in the downward flow region and a minimum in the upward flow one. In the work of Druzgalski et al., [START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] the authors studied, in the chaotic regime, the fluctuating transversal flux averaged in time and in the z direction. Three regions are observed: in the mixing layer, the convection and the electro-migration are dominant; at the overlap of the mixing layer and ESC layer, the electromigration is important; in the quasi-equilibrium EDL and the diffusion layer, the diffusion and the electro-migration are dominant.

In this section and Sec. III F, we consider the two components: their spatial distribution and the local motion of the ionic species in the whole polarization layer. In order to describe the local ionic transfer, the ionic streamlines are computed by means of the ionic velocity defined as U ± = J ± /c ± (Figs. 11 and14). The spatial distribution of the two components of the three relative contributions is displayed in Figs. 12, 13, 15, and 16. The term relative means the ratio of the magnitude of each contribution to the sum of the three ones for the two components

Em ± x = J ± Em x J ± Fi x + J ± Em x + J ± Co x
, for example. In the following, the region 0 < z < 0.6 will be used only in order to make the comments easier.

Figure 11 overlays the cationic streamlines with the color plot of the cationic flux amplitude for the two components. The In the region x > 0.6, the cationic rolls form a barrier for the cationic transfer to the membrane so that the flux is displaced to the outer boundary of the domain and they push the transversal concentration gradient to the region x > 0.3 [Fig. 12(c)]. This displacement induces an increase of the gradient amplitude. Therefore, around the symmetry axis

J + Fi x
is dominant, whereas the region where J + Em x controls the cationic transport is moved to the outer boundary (z = 0 and z = 1.2) [Fig. 12(b)]. The upward displacement of the concentration gradient means that the rolls separate the ESC layer from the diffusion layer and control the maximum expansion of the ESC layer.

In the region 0.3 < x < 0.6, J + = 0 (0.2 < z < 0.4), the transversal motion is controlled by the electric field. Around the symmetry axis, the cationic convection is in the upward direction ( J + Co z < 0) and is opposed to the diffusion ( J + Fi x > 0) and to the electro-migration ( J + Em x > 0). However, the convection is not the dominant mechanism.

In the region 0.1 < x < 0. where J + Fi x cannot be neglected. This layer is located at the transition between the diffusion layer and the ESC one. In these two layers, the concentration gradient is small but between them a jump of concentration takes place over a small distance.

In the vortex region, the negative contribution of J +

F. Description of the anionic transport

Figure 14 displays anionic vortices centered at (x = 0.3, z = 0.2, and z = 1.). Figures 15 and16 show that the spatial distribution of the relative anionic flux contributions is almost the same as the cationic one. This is surprising if we consider the difference between the cationic and the anionic streamlines. But this time, the mean flux is null over the computational domain, J -

G. Potential and pressure disequilibrium on both sides of a vortex

In Sec. III C, we have observed that the fluid motion is induced by the electric and the pressure forces distributed on several layers inside the non-equilibrium EDL and the diffusion layer. We are now able to describe how the distribution of the charge density induces a macroscopic drop of pressure and of potential between z = z 1 = 0 and z = z 2 = 0.6. The idea commented briefly in the work of Pham et al. [START_REF] Pham | Direct numerical simulation of electroconvective instability and hysteretic current-voltage response of a permselective membrane[END_REF] is thereby used in this section but developed in detail. Thus the longitudinal disequilibrium of pressure and potential, leading to the electrokinetic instability, is analyzed quantitatively and compared to the spatial distributions of F Pr z and of F E z . We also link the longitudinal pressure drop to the transversal potential profiles.

Let us describe first the charge distribution in the nonequilibrium EDL at z 1 and z 2 . The non-equilibrium EDL has a smaller width at z 1 than at z 2 and the charge density is much higher at z 1 than at z 2 in the ESC layer. At z 1 , the charge density profiles have a minimum value Q min = 0.016 at x min = 0.0136, a maximum value Q max = 0.0217 at x max = 0.0342, and the position where Q Q max = 1/5 is x = 0.06. At z 2 , Q min = 0.0014 at x min = 0.05, Q max = 0.0015 at x max = 0.12, and Q/Q max = 1/5 at x = 0.26. The value of the total electric charge contained over all the transversal direction ( 1 ∫ 0 Q (x, z) dx) varies between 0.22 and 0.2 as z increases from z 1 to z 2 . At z = z 1 , the total electric charge in the ESC layer is half the electric charge contained in the quasi-equilibrium EDL. At z 2 , the value of the charge ratio falls to 0.07. Therefore, owing to the invariance of the total charge, the width of the quasi-equilibrium EDL, i.e., the charge density, at z 2 must be greater than the quasiequilibrium EDL width at z 1 .

Let us consider now the potential drop between the two axes z 1 and z 2 . We have observed that ∇Φ z << ∇Φ x in the non-equilibrium EDL. So, the Poisson equation becomes -∂ 2 ∂x 2 Φ ∼ Q. By integration, the equation has the following form: (21) where E is the electric field. This means that the variation of the potential gradient or in other word the slope variation of the potential profile between two points is equal to the electric charge contained in this interval. Figure 17(a) displays the potential profiles at z 1 and z 2 . Three regions can be observed. In the first one, located at 0 < x < 0.005, i.e., inside the quasi-equilibrium EDL, the two profiles are very close because Φ = 0 at the membrane surface and Q (0,0.005) is equal to 0.0012 and 0.0014 at z 1 and z 2 , respectively. The two curves are almost tangential at the membrane interface. As x increases from 0 to 0.005, the gradient decreases from 1100 to 500 and from 980 to 300 at z 1 and z 2 , respectively. So the gradient difference is slightly higher at z 2 than at z 1 , but the ratio of the gradient at z 1 to the gradient at z 2 is greater at x = 0.005 than at x = 0. However, over a small distance it has a small effect on the potential profiles.

E x (z, x 1 ) -E x (z, x 0 ) ∼ x 1 x 0 Q (z, x)dx = Q (x 0 ,x 1 ) (z) ,
A 2nd region with a smaller slope is located between 0.005 and x with x = 0.04 and 0.125 at z 1 and at z 2 , respectively. These transversal positions x are slightly higher than the position of Q max . Q (0.005,0.04) (z 1 ) = 7.1 × 10 4 > Q (0.005,0.125) (z 2 ) = 3.6 × 10 4 contrary to the 1st region. Therefore, as x increases over the region II, the difference in gradient is greater at z 1 than at z 2 , i.e., it decreases from 500 to 150 and from 300 to 82, respectively. But Q (0.005,0.04) (z 1 ) is not high enough for the gradient to decrease sufficiently and for the two curves to intersect. Even more, along the 2nd region at z 1 , the potential difference between the two longitudinal positions increases and becomes noticeable. A 3rd region is located between x and 1 with x = 0.08 and 0.2 at z 1 and z 2 , respectively. It is characterized by a much smaller slope than in the 2nd region. We must recall that the term "much smaller" has to be understood in the relative manner. The gradient at x is equal to 26.4 and 18 at z 1 and z 2 , respectively. This 3rd region is located outside the ESC layer within Rubinstein's meaning, i.e., x ≥ x o . In the transition zone (2/3) located between x max and the outer edge of the ESC layer, Q

(x ,x ) (z 1 ) = 2.3 × 10 -4 > Q (x ,x ) (z 2 ) = 1.3 × 10 -4
. So combining the fact that the charge Q (x ,x ) (z 1 ) is high enough to decrease the potential gradient in a significant amount and the fact that x is much smaller at z 1 than at z 2 , the two curves intersect at x = 0.15. This position corresponds to the mean position of the outer limit of the 3rd layer, mentioned in Sec. III C about the transversal gradients, where also F Pr z < 0 and F E z > 0. So, the region (x > 0.15) is characterized by Φ(z 2 ) > Φ(z 1 ). This coincides to the presence of a negative longitudinal electric force around the symmetry axis, inducing with the pressure force the divergence of the fluid flow [also see Figs. 13(b) and 16(b)]. The maximum of potential difference is located at x = 0.04. This corresponds to the mean position of the maximum of F E z . Figure 17(b) displays the transversal pressure profiles at z 1 and z 2 . The inset of Fig. 17(b) is a magnification of the profiles in the non-equilibrium EDL. We have observed in Sec. III D that the transversal mechanical equilibrium is fulfilled in this layer. Therefore, from Eq. ( 20), P(x, z) ∼ Pe/2 E 2 x . Contrary to the potential, the pressure at the membrane surface is extrapolated. Its value is equal to 2.4 × 10 5 and 1.95 × 10 5 at z 1 and z 2 , respectively. P(z 1 ) is greater than P(z 2 ) but the relative difference is small because the transversal potential gradient is slightly higher at z 1 than at z 2 at the surface membrane.

The 1st region (x < 0.005) is located in the quasiequilibrium EDL. It is characterized by a high decrease in the pressure because in this region ∇Φ x , i.e., the electric field magnitude, decreases a lot as discussed above. For the same reason about the difference in potential gradient between the two longitudinal positions, we observe P(z 1 ) > P(z 2 ).

The 2nd region is characterized by a smaller pressure gradient which is almost constant. This region lies between 0.01 and 0.035 and between 0.03 and 0.1 at z 1 and z 2 , respectively. In Fig. 17(a), this pressure region is located just before the slope change at the transition region 2/3. In the ESC layer, it takes place between x min and x max . So, considering the mechanical equilibrium, the pressure gradient is constant because the increase in charge density balances the decrease in potential gradient. Figure 17(a) shows that in the transition between regions 2 and 3 at z 1 , the potential gradient decreases noticeably and E x (z 1 ) < E x (z 2 ) if x > 0.04. This is why the pressure curves intersect at x = 0.04. From the mechanical equilibrium in the transversal direction, we can infer that ∇P

x (z 1 ) > ∇P x (z 2 ) when x > 0.01 regardless of E x because the electric charge density is very small at z 2 and is the key parameter. This limit at x = 0.04 corresponds to the mean position of the outer boundary of the 2nd layer, mentioned in Sec. III C about the transversal gradients, where also F Pr z > 0 and F E z > 0. Figure 17(b) shows a maximum of the pressure difference at x = 0.06. As x increases from 0.035 to 0.06, the maximum of the longitudinal pressure profile displaces from z = 0.3 (and z = 0.9) to z = 0.5 (and z = 0.7). In the region 0.06 < x < 0.2, the maximum of the longitudinal profile is located at z = 0.6. The pressure plateau, mentioned in Sec. III C, appears when x > 0.2.

The 3rd region, characterized by a small pressure decrease, is located beyond the outer edge of the ESC layer: x > 0.08 at z 1 and x > 0.3 at z 2 . The pressure at z 2 remains greater than the pressure at z 1 . From Fig. 17(a), we can deduce that this is the result of the compression of the ESC layer at z 1 and the result of a smaller potential gradient at z 1 than at z 2 in the transition region 2/3.

IV. CONCLUSIONS

In this work, ion transfer and hydrodynamics are studied in the presence of a steady-state electro-kinetic instability at the surface of a cation-exchange membrane. In this aim, the Poisson-Nernst-Planck equations coupled with the Stokes equations are solved in the polarization layer until the Debye scale. After a validation of the numerical method, the results are presented in three parts.

In the first part, marginal instability properties are studied by varying the cation concentration at the fluid/membrane interface c + interf and by investigating the stability of the vortex pair size as a function of the imposed potential drop ∆Φ. The numerical results show that the marginal instability threshold ∆Φ * increases with c + interf and the current density I is a decreasing function of c + interf . ∆Φ * and I are logarithmic functions of c + interf . This is in opposition with the theoretical results published by Rubinstein and Zaltzman. [START_REF] Rubinstein | Electro-osmotic slip of the second kind and instability in concentration polarization at electrodialysis membranes[END_REF][START_REF] Zaltzman | Electro-osmotic slip and electro-convective instability[END_REF] The asymptotic development leads to an expression of the slip velocity independent on c + interf . So, the simplified model describing the electro-kinetic instability in the diffusion layer cannot predict the dependency of the instability on c + interf . A higher order development may be necessary to confirm the numerical results. However, the high difference in velocity value between the analytical expression and the numerical computation remains unexplained. As regards the vortex pair size, the numerical results carried out with L z = 6 and c + interf = 2 show that close to the instability threshold the pair size value of L z /5 is the most stable in agreement with the published results. [START_REF] Demekhin | Direct numerical simulation of electrokinetic instability and transition to chaos motion[END_REF][START_REF] Druzgalski | Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface[END_REF] However, the size L z /4 is also stable. As ∆Φ increases, the range of stable vortex pair size increases so that at ∆Φ = 30, L z /6 and L z /7 seem to be the most stable sizes.

The second part deals with the spatial distribution of the hydrodynamic and of the ion transfer. As concerns the hydrodynamic, the fluid motion along the membrane is mainly induced by the electric force located around the maximum value of the electric charge density in the ESC region all along the membrane. But inside the quasi-equilibrium EDL and around the inner edge of the ESC layer, the fluid motion is also induced by the pressure gradient so that the fluid velocity increases as the fluid approaches the symmetry axis showing a kinetic energy spot in the region where the ESC layer width increases. Owing to the incompressibility, the upward flow displays a spot of kinetic energy in the region between the counter-rotating vortices where the streamlines are no more convergent. In this region, the pressure force, directed upward, is far greater than the electric one. The pressure does not depend on z contrary to U x which is a quadratic function of z. The divergence of the velocity field is mainly induced by the longitudinal pressure and electric forces directed to the outward direction.

In the analysis of ion transfer, three contributions to the ionic velocity (ratio of the ionic flux to the ionic concentration) are compared: the Fickian diffusion, the electro-migration, and the convection. As concerns the cations, the transfer is mainly directed to the membrane surface due to the electric force and the convection successively as we move away from the boundary z = 0. However, a stagnation zone takes place between the counter-rotating hydrodynamic vortices, i.e., where the flow is upward. The small ionic velocity is the result of the transversal balance between the convective flux and the electro-migration one. In the two counter-rotating cationic vortices, the longitudinal motion is induced by the Fickian diffusion at the bottom of the vortices and by the electro-migration at the top. As concerns the anions, the numerical results show two regions between which there is no anionic flux. The first one is located against the stationary reservoir. The anions enter the polarization layer between the hydrodynamic vortices by diffusion and are pushed out of the layer in the outer boundary region by the electric force. In the second region, lying between the membrane and the first zone, the anionic motion is characterized by two counter-rotating vortices. The longitudinal transfer is induced by convection and Fickian diffusion. In the transversal direction, the electro-migration and the convection induce the downward motion one and the electro-migration controls the upward one.

Finally, using the knowledge achieved from the spatial analysis of the fluid motion and of the ionic fluxes, we describe, in another way, how the spatial charge distribution, in the quasiequilibrium EDL and in the ESC layer, and also the size of these two layers influence the longitudinal potential and pressure difference. In particular, we are able to explain at the macroscopic scale the presence of two layers along the membrane surface where the pressure force and the electric force induce the fluid motion underlying the electro-kinetic instability.

FIG. 1 .

 1 FIG. 1. Schematic representation of the numerical model. Domain dimension: L x = 1 and L z = 6. The boundary conditions are indicated in brackets. The color plot represents the electric charge density. The streamlines visualize the instability. The red arrows indicate the direction of the vortices rotation. ν = 10 3 , ∆Φ = 20, c + interf = 2, c o = 1.

2 z

 2 Ω converges to the value of 10 13 after an iteration number (N iter ) of 4 × 10 5 . Using this new solution for ∆Φ = 19.35, U 2 z Ω increases exponentially with time if the iteration number is greater than 3.5 × 10 5 . At N iter = 9 × 10 5 , U 2 z Ω reaches the value of 10 7 . So, the potential drop threshold in the increasing direction lies between 19.25 and 19.

FIG. 5 .

 5 FIG. 5. Current density versus electric potential drop for two values of domain length L z and several values of c + interf . Solid lines: computations carried out in the case of a quiescent solution.

FIG. 6 .

 6 FIG. 6. Effect of c + interf on: (a) the electric potential threshold ∆Φ * , (b): the current density. Two values of domain length L z : 6/5 and 6/7. ν = 10 3 . Symbols: numerical results. Solid lines: analytical expression.

  density (F = I) for two electric potential drops and two domain lengths. ν = 10 3 . L z a (∆Φ = 22.5) b (∆Φ = 22.5) a (∆Φ = 25) b (∆Φ = 25)

.

  FIG. 7. Variation of the amplitude of several wavelengths vs time. ν = 10 3 , c + interf = 2, and L z = 6. (a): ∆Φ = 20, (b): ∆Φ = 30.

  order to investigate deeper the stability of the vortex pair as a function of its size, numerical experiments are carried out in two steps with the same parameter and ∆Φ values. This time, the stability is studied in the non-linear domain by means of a high amplitude of perturbation in the 2nd step. The 1st step consists in computing the numerical solution with a 1D initial solution at ∆Φ = 20 and 30 to which is added a small amplitude monochromatic sinusoidal signal in the z direction (λ o = L z /N o = L z /4, L z 5, L z /6, L z /7). This allows us to obtain a vortex pair size of λ o . In the 2nd step, a high amplitude monochromatic sinusoidal signal is added to the vortex array. The wavelength (noted as λ o = L z /N o ) is different from λ o and the amplitude is this time 20%. The computation is stopped when the anionic flux converges to zero. The vortex pair size reaches the final value λ final . The time evolution of the roll structure is not taken into account. At ∆Φ = 20 the 1st step, the steady instability is reached with a vortex pair size equal to the signal wavelength λ o but the length ratio N o = L z /λ o = 7. With this initial value, the stabilized vortex pair size is L z /4. For the 2nd step, the results in Table II show that the vortex pair sizes L z 5 and L z /4 are the most stable. It must be added if λ o = L z /6 and λ o = L z /5, the final vortex pair size reaches the value of L z /5 after the intermediate value of L z /4.At ∆Φ = 30 the 1st step, the vortex pair size is different from λ o when λ o = L z 4 only. In this case, the vortex pair size is L z /6. For the 2nd step with λ o = L z /5, L z 6, L z /7, if the perturbation has a wavelength smaller than the initial vortex pair size (λ o = L z /N o , λ o = L z /(N o + 1)), the final size moves to λ o except λ o = L z /7. On contrary, if the perturbation wavelength is greater than the initial vortex pair size (λ o = L z /N o , λ o = L z /(N o -1)), the vortex pair size remains unchanged except λ o = L z /8. This confirms that the stability of a vortex pair decreases with its size.

  FIG. 8. Color plot of the electric charge density Q = c + c (a) and of the kinetic energy (b) in the polarization layer. The counter-rotating vortices are visualized by the streamlines. ∆Φ = 25, ν = 10 3 , c + interf = 1, L z = 6/5. The red arrows indicate the direction of the vortices rotation.

FIG. 9 .

 9 FIG. 9. Color plot of the difference in magnitude between the electric force and the pressure force. (a): Longitudinal component. (b): Transversal component. The counter-rotating vortices are visualized by the streamlines. ∆Φ = 25, ν = 10 3 , c + interf = 1, L z = 6/5. The red arrows indicate the direction of the vortices rotation.

  gathers the values of x o , the estimated position x o of the minimum of electric charge density at the edge between the ESC region and the quasi-equilibrium EDL, and the ratio of the numerical velocity to the theoretical slip velocity. These values are computed for three longitudinal positions where E c reaches the highest value. The parameter values used in the computation are ∆Φ = 25, c + interf = 1, ν = 10 3 , and L z = 6/5. The values of x o are in agreement with the assumption of outer edge of the ESC layer. At these transversal positions, the value of c + is around 1/5 of the maximum value located between 0.04 and 0.066. The position of the minimum of concentration x c + min is exactly estimated by the analytical relation found by Urtenov et al.,[START_REF] Urtenov | Basic mathematic model of overlimiting transfer enhanced by electroconvection in flow-through electrodialysis membrane cells[END_REF]

Figure 10

 10 FIG. 10. Transversal profiles for three longitudinal positions. Open circle: plot of the longitudinal velocity component. Solid lines: plot of the 2nd order polynomial fit. ∆Φ = 25, ν = 10 3 , c + interf = 1, L z = 6/5.

  last column). At x = x o , the estimated slip velocity is equal to 586, 1100, and 650 as z increases from 0.3 to 0.5. So these values are far much greater than the maximum value of the numerical velocity and do not have the same sign at z = 0.5. Comparison with smaller values of ν until 10 4 , higher values of c + interf until 4 and at ∆Φ = 20 also leads to the same disagreement.

FIG. 11 .

 11 FIG. 11. Color plot of the transversal component of the cationic flux (a) and of the longitudinal one (b). The streamlines display the cationic motion. ∆Φ = 25, ν = 10 3 , c + interf = 1, L z = 6/5. The blue arrows indicate the direction of the vortices rotation.

  cationic motion, and the cations move to the outer domain boundary by convection. Owing to the fluid vortices, this region may be divided in three parts in the z direction [Figs. 12(a) and 12(b)]. Close to the outer boundary (z < 0.2), the cations move to the membrane by convection. Around a curve joining the points where U x

  FIG. 13. Color plot of the relative amplitude of each cationic flux contribution. Longitudinal component. (a): Convective flux, (b): electro-migration flux, (c): diffusion flux. The streamlines represent the fluid motion. ∆Φ = 25, ν = 10 3 , c + interf = 1, L z = 6/5. The red arrows indicate the direction of the vortices rotation.

  . Added to the J + Co z one, their contribution is close but greater in magnitude than the J + Fi z one if x > 0.25, i.e., in the upper part of the vortices, and push the cations in the outward direction. If x < 0.25 (lower part of the vortices), J + Fi z is again dominant and the cationic streamlines are convergent.

  FIG. 15. Color plot of the relative amplitude of each anionic flux contribution. Transversal component. (a): Convective flux, (b): electro-migration flux, (c): diffusion flux. The streamlines represent the fluid motion. ∆Φ = 25, ν = 10 3 , c + interf = 1, L z = 6/5. The red arrows indicate the direction of the vortices rotation.

FIG. 16 .

 16 FIG. 16. Color plot of the relative amplitude of each anionic flux contribution. Longitudinal component. (a): Convective flux, (b): electro-migration flux, (c): diffusion flux. The streamlines represent the fluid motion. ∆Φ = 25, ν = 10 3 , c + interf = 1, L z = 6/5. The red arrows indicate the direction of the vortices rotation.

  (a), and16(c) show that the anions move mainly by diffusion in the longitudinal direction when 0 < x < 0.25.

FIG. 17 .

 17 FIG. 17. Transversal profiles for two longitudinal positions. (a) Plot of the electric potential ∆Φ(x), (b): plot of the pressure. ∆Φ = 25, ν = 10 3 , c + interf = 1, L z = 6/5.

  At the same time, the concentration gradient at the outer edge of the quasi-equilibrium EDL appears. In a second step (t ∈ [0.01, 0.04] or ν I (t)/2 ∈ [0.4, 0.063]), the ESC region appears and the diffusion layer spreads to the bulk boundary (region II). During the third one (t ∈ [0.04, 0.2] or ν I (t)/2 ∈ [0.063, 0.016]), the ESC region remains unchanged, the diffusion layer continues to spread,

	and U	2 z Ω

16 is close to I lim as noticed in Ref. 29. Four regimes can be observed for the four values of ∆Φ. For example, when ∆Φ = 30, during the first step (t ∈ [0, 0.01] or ν I (t)/2 ∈ [5, 0.4]), the quasi-equilibrium EDL FIG. 4. Temporal variation of the current density for several values of ∆Φ (solid lines). Symbol: self-similar curve. Four regimes are displayed for ∆Φ = 30. ν = 10 3 , c + interf = 2.

develops (region I).

  d is negative in the case of a cation-exchange membrane. Its magnitude increases with the fixed charge density. Assuming that c + interf ∼ Y /2, as usually observed at equilibrium, Φ d reaches a value of 2.77 and 4.16 if c + interf = 8 and 32, respectively, and the interfacial potential has the value of 2.08 and 3.47

TABLE II .

 II Dependency of the final vortex pair size (λ final ) as a function of the wavelength perturbation (λ o ). λ o : Initial vortex size. ∆Φ = 20, ν = 10 3 , c + interf = 2, L z = 6.

  ) (U z (x o )/U s (x o )) at x o .

	z	x o	x c + min	x o	U z (x o )/U s (x o )
	0.3	0.9	0.90	0.07	0.053
	0.4	0.979	0.98	0.08	0.024
	0.5	0.97	0.97	0.14	0.001 51


Wang Fan

Samia Ainouz

Fabrice Meriaudeau

Abdelaziz Bensrhair

Abdelaziz Bensrhair Polarization

POLARIZATION-BASED CAR DETECTION

Keywords: Car detection, polarization, feature selection, DPM, road scenes

Road scene understanding is a vital task for driving assistance systems. Robust vehicle detection is a precondition for diverse applications particularly for obstacle avoidance and secure navigation. Color images provide limited information about the physical properties of the object. This results in unstable vehicle detection caused mainly from road scene complexity (strong reflexions, noises and radiometric distortions). Instead, polarimetric images, characteristic of the light wave, can robustly describe important physical properties of the object (e.g., the surface geometric structure, material and roughness etc). This modality gives rich physical informations which could be complementary to classical color images features. In order to improve the robustness of the vehicle detection purpose, we propose in this paper a fusion model using polarization information and color image attributes. Our method is based on a feature selection procedure to get the most informative polarization feature and color-based ones. The proposed method, based on the Deformable Part based Models (DPM), has been evaluated on our self-collected database, showing good performances and encouraging results about the use of the polarimetric modality for road scenes analysis.

INTRODUCTION

For an intelligent device in road scenes, cars appear to be one of the most frequently observed yet dangerous objects. Car detection has broad range of applications such as autonomous driving or obstacle detection and avoidance. It is a challenging problem due to its large structural and appearance variations.

Dalal et al. [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] proposed a single filter based on Histogram of Oriented Gradients (HOG) features to represent an object. The model is a single filter that slides throughout the image. A score is computed with each position and scale. Based on this effective feature, Felzenswalb et al. [START_REF] Felzenszwalb | Object detection with discriminatively trained part-based models[END_REF] proposed Deformable Part Models (DPM). The object is modeled by several deformable parts to better handle the object appearance points of view (frontal, rear, left side or right side) and to inte-grate the intra-class variations. Specifically for car detection, Wu [START_REF] Wu | Learning andor models to represent context and occlusion for car detection and viewpoint estimation[END_REF] proposed a reconfigurable hierarchical and/or model to integrate the context and the occlusion patterns, in which the DPM is utilized as the deformable feature.

Classical color-based detection methods [START_REF] Sun | Onroad vehicle detection: A review[END_REF], [START_REF] Mukhtar | Vehicle detection techniques for collision avoidance systems: A review[END_REF], [START_REF] Garcia | Sensor fusion methodology for vehicle detection[END_REF] including those mentioned above, depend extremely on outdoor illumination conditions. These conditions, in the road scene context, are most of the time hard and complex (strong reflexions, presence of transparent objects, occlusions, bad weather conditions) making RGB images poor and limited [START_REF] Wang | Polarization-based specularity removal method with global energy minimization[END_REF].

The polarization of the reflected light is strongly linked to the physical properties of the surface [START_REF] Wolff | Polarization-based material classification from specular reflection[END_REF]. Compared to other sensors, polarimetric imaging systems have the advantage to characterize the geometric aspect of the surface (orientation, reflexion angle, degree of polarization), to identify the material nature of each object (refraction index, material type) and to offer physical and optical information (reflexion, transmittance, depolarisance). Polarization is able to bring additional information of an object other than its intensity [START_REF] Lawrence | Surface orientation from polarization images[END_REF].

In this paper, we propose to use polarization features as a complementary information to color-based ones in order to improve car detection results. To our knowledge, this is the first work in the literature that attempts to use polarizationbased features for outdoor object detection. A feature selection process is performed to select the most informative polarization feature. The detection scheme is carried out using DPM detector. Polarization-based DPM detector and a color-based DPM detector are trained independently and different score maps are produced by the two models. A fusion rule that takes the polar-based model as a confirmation to the color-based one is proposed to achieve the final detection. Experiments performed on our self collected dataset, show that fusing polarization and color features reduces strongly false alarm rate (false bounding boxes), and improve effectively the car detection accuracy.

POLARIZATION FORMALISM

The main applications of polarization imaging are related to distinguish metallic object from dielectric surface [START_REF] Born | Principles of optics: electromagnetic theory of propagation, interference and diffraction of light[END_REF], [START_REF] Bass | Handbook of optics[END_REF]. Polarization imaging permits likewise to give threedimensional information of specular [START_REF] Morel | Calibration of catadioptric sensors by polarization images[END_REF], [START_REF] Wolff | Constraining objects feature using a polarization reflectance model[END_REF] and transparent objects [START_REF] Ferraton | Optimization of a polarization imaging system for 3d measurements of transparent objects[END_REF]. The physical principle of the polarization of the light is that after being reflected, an unpolarized light wave become partially linearly polarized depending on the surface normal and on the refractive index of the material it impinges on [START_REF] Wolff | Polarization-based material classification from specular reflection[END_REF], [START_REF] Lawrence | Surface orientation from polarization images[END_REF], [START_REF] Morel | Polarization imaging applied to 3d recontsruction of specular metallic surfaces[END_REF] . The reflected partially linearly polarized light is described by a measurable vector, named the Stokes vector S = [S 0 , S 1 , S 2 ]. The first component S 0 relates to the object intensity, the two others describe the physical properties of the object. From these three components, other physical properties are derived such as the light magnitude I , the degree of polarization (DOP) ⇢ and the angle of polarization (AOP) ' [START_REF] Ainouz | Adaptive processing of catadioptric images using polarization imaging: towards a pola-catadioptric model[END_REF]. In order to measure the polarization parameters, at least three images are required. For this purpose, a rotating linear polarizer around three angles (↵ i ) i=1:3 is placed in front of the camera. Figure .1 gives an example of a polarization device. For each angle (↵ i ), an intensity I(↵ i ) is measured by the camera. The relationship between the acquired images and the Stokes vectors is given by :

I(↵ i ) = 1 2 [1, cos(2↵ i ), sin(2↵ i )].[S 0 , S 1 , S 2 ] t (1)
The DOP is calculated as

⇢ = p S 2 1 +S 2 2 S0
and the AOP as

' = 1 2 tan 1 S2 S1 3. FEATURE EXTRACTION 3.1. Polarization feature selection
Because of the noisy nature of the polarization parameters, a feature selection is required to find the most informative polar-based feature for the detection model training. In our case, the feature selection procedure is trained by the Dalal-Trigg detector [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] by replacing the HOG feature by the polarization features.

To train the Dalal-Trigg detector, similar to the original HOG features, the polarization features are extracted based on blocks and cells. A 8 ⇥ 8 block was divided into four 4 ⇥ 4 cells. For each pixel in the cell, a feature vector is extracted, which contains the 3-dimensional Stokes vector, the DOP (⇢) and the AOP ('). The feature vector [S 0 , S 1 , S 2 , ⇢, '] of a By applying the feature selection process presented above, the AOP was selected as the most informative feature with regards to the car detection purpose.

The Angle of polarization

The AOP refers to the direction of the polarization of the reflected light. It is determined by the angle of the incident light (generally for outdoor applications, the incident light is assumed to be unpolarized), the surface orientation of the object and the material of the object. For rough surface, as surface orientations of neighboring pixels change a lot, the AOP changes in an irregularly way. For smooth surfaces, however, the AOP changes smoothly and continuously. Up to the road scenes analysis, especially for car detection tasks, the AOP on the car surfaces changes gradually according to the surface geometry structure, while it is noisy for other objects. This phenomena is illustrated in Figure . 2 where the AOP image on the tree area is highly noisy, on the road it is better but still much noisy than on the car. It can be observed that the AOP image describes the geometry structure of the car, which is even more clear than that from the color image.

We can conclude that, the AOP describes the geometry structure that cannot be captured from the color image. In addition, it highlights the area of the car among the noise.

DETECTION MODEL

The DPM model proposed by Felzenswalb [START_REF] Felzenszwalb | Object detection with discriminatively trained part-based models[END_REF] is known to be an efficient method that handles the intraclass variations of an object. The DPM is a filter-based detector applied throughout a feature map to get a detection score. A feature map is an array that each entry is a d dimensional feature vector computed from the corresponding image location. Let be the feature map, with (x, y) corresponds to the entry at the location (x, y). The filter f shares the same dimension as the feature map, and the detection score is computed as the dot product f • (x, y). The detection score for a given location is proportional to the possibility of the existence of an object on that location (refer to Felzenswalb [START_REF] Felzenszwalb | Object detection with discriminatively trained part-based models[END_REF] for more details about the initialization and optimization process).

FUSION RULE

As polarization provides complementary information that is not accessible from color images, the fusion of polarization and color improves the result provided by color-only methods [START_REF] Martin | Polarization distance: a framework for modeling object detection by polarization vision systems[END_REF] [START_REF] Lawrence | Polarization vision: a new sensory approach to image understanding[END_REF]. In order to improve the detection results, we have proposed a fusion rule by taking the polarization result as a confirmation to the color one.

Following the training process, we obtained a polardeformable part model M p using the AOP image and a colordeformable part model M c using the color image. Since each location is scored by the model, a score map is formulated for each domain (polarimetric and color). The higher score a location gets, the strong the possibility that a car appears at this location. The result is given by simply thresholding the score map by a trained threshold. The idea is that the score location on which M c produces high scores are confirmed by the score of M p . If both M c and M p produces high score, this location is then considered as an object, and vice versa.

The score map given respectively by M c and M p are shown with the original image in Figure . 3. It can be seen that on the color-based score map, both the car area and the road area get high scores, which consequently generates false positive.

The score M p is used as a confirmation to M p in order

Results

As a standard evaluation, the detected bounding boxes on the test set was assessed by the evaluation scheme provided by the PASCAL VOC data set [START_REF] Everingham | The pascal visual object classes (voc) challenge[END_REF].

The Precision-Recall Curve (PRC) is firstly computed and shown in Figure d + color (fusion). The corresponding Average Precision (AP) of each curve is shown in Table 1.

It can be observed from both Figure. [START_REF] Sun | Onroad vehicle detection: A review[END_REF] and Table 1 that the fusion of AOP with color produces the best result, followed by the fusion of 20 d + color which performs only slightly better than the color only one. The 20 d or the AOP alone This result confirms that once the AOP is properly fused with the color model, it provides complementary information that improves the AP by 3.4%. The 20 d is more stable than the AOP, however, with all the redundant features and the more complex model, it almost does not improve the colorbased result (0.1% can be even neglected). Moreover, the AOP alone achieves 31.8% of AP, which is important comparing to the 53.8% performed by the whole polarization feature 20 d. This result shows that the AOP has an important impact on the detection process. According to the above analysis, it can be concluded that our proposed pre-selection method is valid. By using the selected feature (AOP) and the proposed fusion rule, the polarization features provided useful information which improved and reinforced the color-based method.

To evaluate the improvement provided by the fusion of AOP + color, we compare the detected bounding box by AOP + color method with the color-only one. This comparison is shown in Figure . 5, The first and third rows show the results of the detection with the AOP + color fusion and the second and fourth rows show the corresponding color-only detection method. It is worth to note that, after fusing the two different sources of information via the proposed fusion scheme, while keeping the true positive bounding boxes, the false detection are effectively removed.

DISCUSSION

The detection results have shown that the proposed methodology largely reduced the false detection rate and enhance the robustness of the model. Our approach confirmed that polarization-based features can provide useful information for car detection. False alarms were reduced by a simple but effective fusion rule. However, as it can be seen in Figure . 6, the polarization-based feature model does not generate new true positive detection boxes. This limitation is caused by the And-fusion used scheme, which is too much selective.

Our first objective by this work is to prove that polarization is the suitable alternative to color-based models for improving the detection result. This first observation is encouraging to continue our investigation in this direction. As a future work, the polarization feature should be properly fused with the color one inside a training loop. A stable model which integrates the color-based feature and the polar-based feature by an early-based fusion scheme might be able to both reduce the false detection and produce new true positive bounding boxes as well. Other more recent methods like the transfer learning based on the Deep learning networks, that shown interesting performances in several computer vision applications even for vehicle detection [START_REF] Chabot | Deep manta: A coarse to fine many task network for joint 2d and 3d vehicle analysis from monocular image[END_REF], should be tested in the polarization domain.

The enhancing of the polarization images quality is not to be missed. Because of the noise reaching polarization aquisitions, it was difficult to leverage from the whole polarimetric information. For example, the DOP coupled with the AOP could give better results in the detection process provided that the DOP is properly calculated. An adapted (physical) filtering process is thus necessary to get better results.

The last point that deserves a discussion in this work is the real time constraint. The proposed algorithm takes around 5 seconds for the detection task on an image of 320 ⇥ 240, without tacking into account the learning step time. It is far away to be a real time achievement. The use of speed up tools like Open CV implementation and GPU configuration should make faster the detection task .

Fig. 1 .

 1 Fig. 1. The polarization device principle.

Fig. 2 .

 2 Fig. 2. The original image and its corresponding AOP.

Fig. 3 .

 3 Fig. 3. The score map given by color-based and AOP-based methods.

Fig. 4 .

 4 Fig. 4. Precision recall curve on the test set

1 .

 1 Data-set and configurationSince there does not yet exist a public polarization-based data set, our experiment was implemented on a self-collected data set1 . It contains 153 scenes, with three images I 0 , I 45 and I 90 for each scene taken manually by the corresponding polarizer angles. All the cars in the image were labeled by bounding boxes. We divided the data set so that 115 scenes were used for the training and 38 scenes for the model evaluation. The parameters of the DPM model are set as proposed by Felzenszwalb et al[START_REF] Felzenszwalb | Object detection with discriminatively trained part-based models[END_REF].

  .

4 .

 4 The PRC of respectively the AOPbased model, the color-based model and the fusion of AOP with color are shown in the same figure. We trained a model with the 20-dimensional polarization features (the features described in section 3), fused with the color based model. The results are noted in Figure. 4 by 20 d (polar-based model) and 20

Fig. 5 .

 5 Fig. 5. Detection results. The first and third rows refer to the polar-based model whereas the second and fourth rows are the corresponding color-only based model.

Fig. 6 .

 6 Fig. 6. Main limitations. False alarms are reduced however not new true positive detection boxes are generated.

Table 1 .

 1 The Average Precision (AP) rate does not provide comparable results, this is because the noise still strongly present in polarization acquisitions.

http://pagesperso.litislab.fr/fwang/fichiers/