Jean-Pierre Simonin
email: jpsimonin@gmail.com

José Miguel Ramos

José Torres-Arenas

Diffusion coupling in multiply associating electrolyte solutions

Keywords: Electrolytes, association, non-ideal, diffusion, mean-spherical approximation (MSA)

Diffusion coupling in multiply associating aqueous electrolyte solution is studied. The process is illustrated in the case of zinc(II) ion which forms a set of complexes, ZnCl 2-n n (with n = 1, .., 4), in the presence of chloride ions. The transport of zinc ion, taken in radioactive form 65 Zn 2+ , in the electric field created by a gradient of LiCl, is investigated experimentally by employing an adaptation of the closed capillary technique. The transient diffusion process is modeled by using two different treatments: finite difference (FD) simulation and normal mode (NM) analytic solution. Deviations from ideality are taken into account with the use of the mean spherical approximation (MSA), and the internal electric field is calculated using the dynamical electroneutrality condition. The theoretical FD and NM results are compared with the experimental data about the diffusion of zinc.

Introduction

Aqueous electrolytes are ubiquitous on earth. They may be found in natural waters (e.g., in oceans, lakes, sediments), in most living beings and in plants, trees,... They have a huge influence on the development and functioning of life on the planet. The motion of ionic species governs the local amount of ions. In this respect, aqueous electrolytes have the particular property that the movements of the ions in a medium are strongly coupled through the effect of electrostatic interactions [START_REF] Tyrrell | Diffusion in Liquids[END_REF][START_REF] Mills | Coupling of diffusion Processes in Multicomponent Electrolyte Solutions[END_REF].

Diffusion coupling thus can be observed in many natural and industrial domains and systems, such as in geochemistry (when studying diagenetic fluxes in sediments [START_REF] Yaakov | Diffusion of seawater ions. I. Diffusion of seawater into a dilute solution[END_REF][START_REF] Lasaga | The treatment of multi-component diffusion and ion pairs in diagenetic fluxes[END_REF]), in biophysics (e.g. for a description of permeation through ionic channels of cells [START_REF] Eisenberg | Ionic channels in biological membranes -electrostatic analysis of a natural nanotube[END_REF]), in engineering processes (e.g., when using membranes for industrial separations [START_REF] Simon | Assessment of diffusion coupling effects in membrane separation. Part I. Network thermodynamics modelling[END_REF]), in materials science (e.g., in the study of diffusion of ions in concrete [START_REF] Samson | Calculation of ionic diffusion coefficients on the basis of migration test results[END_REF] or for the assessment of chemical ageing of concrete and composite materials [START_REF] Bai | Advanced Fibre-Reinforced Polymer (FRP) Composites for Structural Applications[END_REF]).

Diffusion coupling in aqueous ionic solutions has been investigated at a fundamental level for a long time [START_REF] Tyrrell | Diffusion in Liquids[END_REF]. Theoretical descriptions have been developed that assume linear transport theory, in which fluxes are linear functions of forces [START_REF] Boudreau | Multicomponent ionic diffusion in porewaters: Coulombic effects revisited[END_REF]. Multicomponent diffusion has been examined in the framework of Onsager formalism of irreversible thermodynamics [START_REF] Miller | Application of Irreversible Thermodynamics to Electrolyte Solutions. I. Determination of Ionic Transport Coefficients l ij for Isothermal Vector Transport Processes in Binary Electrolyte Systems[END_REF].

A phenomenon that has great influence on coupled diffusion in electrolytes is ionic association, consisting either of ion pairing or chemical reaction. In geochemistry, Lasaga first examined the influence of ion pair formation on diagenetic fluxes in marine sediments [START_REF] Lasaga | The treatment of multi-component diffusion and ion pairs in diagenetic fluxes[END_REF]. Ion pairing is important because it modifies the effective charge of the diffusing species and so changes their response to the internal electric field created by the diffusing ions. A clear example is that of a monovalent tracer ion that associates with a monovalent anion and forms a neutral species that is insensitive to the diffusion field [START_REF] Lasaga | The treatment of multi-component diffusion and ion pairs in diagenetic fluxes[END_REF].

Some time ago, we investigated the effect of ion pairing caused by magnesium(II) ion on the transient transport of sulfate ion by carrying out experiments in a special diffusion cell (closed capillary) [START_REF] Simonin | Diffusion coupling in electrolyte solutions. 1. Transient effects on a tracer ion: sulfate[END_REF]. Later we investigated the effect of a pH gradient on sulfate and phosphate ions [START_REF] Simonin | Transport coupling of ions: Influence of ion pairing and pH gradient Application to the study of diagenetic fluxes[END_REF]. In these works, the diffusion-reaction equations were solved by using two types of treatments: a 'normal-mode' analysis in which the equations are linearized and solved analytically [START_REF] Turq | Normal Mode Analysis od Diffusion in Multicomponent Electrolyte Solutions. 1. Nonassociating Electrolytes[END_REF]; and a classical numerical finite-difference method in which time and space are 'sliced' and the equations are solved incrementally in the course of time. In these treatments, the equations were written in the ideal case in which activity coefficients are taken equal to unity. Besides, the associations involved only one type of complex, leading to the ion pair MgSO 0

In environmental media, ionic pollutants exist in various forms. For instance, toxic heavy metal ions of various valencies such as zinc, cadmium, lead or mercury may form a series of complexes with ions like chloride, nitrate, hydroxide,... The speciation of these ions in the environment is an important issue because complexation modifies the electric charge borne by the ions. This phenomenon modifies the transport of trace metal ions when they are submitted to diffusion electric fields caused by major ions in the environmental medium, and also the physical and chemical behavior of the ion interacting with the environment, e.g. their interactions with colloidal and mineral particles.

Various softwares exist commercially which permit estimations of the speciation for many metal ions in aqueous solutions. One may cite the following which have been developed in various countries: MINTEQ [START_REF]Free software[END_REF], MINEQL+ [START_REF]Software developed in the US[END_REF], JCHESS [START_REF]Software developed at Ecoles des Mines in France[END_REF], PHREEQC [START_REF]Free software developed at US Geological Survey[END_REF] and CHEAQS [START_REF] Verweij | Free software developed by[END_REF]. The calculation of the speciation is based on a solution to the chemical association equilibria which involve the introduction of thermodynamic association constants (cumulative constants β n for n = 1, 2, 3, ...) and the use of formulas for the computation of the activity coefficients of the species. The softwares rely on previous determinations of the association constants, which can be found in famous books, e.g., in the book by Sillen and Martell [START_REF] Sillen | Stability Constants of Metal-Ion Complexes[END_REF]) or in the NIST Database [START_REF]NIST Database 46 Version 8.0, NIST Critically Selected Stability Constants of Metal Complexes[END_REF]. In general, deviations from ideality were computed using equations of the Debye-Hückel type, like the Davies equation. Depending on the values taken for the β n 's and on the equation used, softwares may give speciations that differ in magnitude. Because the speciation is generally very difficult to determine experimentally, there is uncertainty in the determination of speciation in electrolyte solutions.

In the present work, we propose to study transient diffusion patterns in the case of a metal cation giving rise to multiple association. The metal cation (in tracer amount) is placed in the electric field created by the diffusion of a major salt. This type of experiment is a variation on the theme of experiments we had carried out with magnesium(II) and sulfate ions [START_REF] Simonin | Diffusion coupling in electrolyte solutions. 1. Transient effects on a tracer ion: sulfate[END_REF][START_REF] Simonin | Transport coupling of ions: Influence of ion pairing and pH gradient Application to the study of diagenetic fluxes[END_REF]. The latter system gave one type of ion pair. In the present work we investigate the case of electrolytes leading to a set of complexes of increasing stoichiometries through a series of stepwise associations.

Moreover, in contrast with our earlier work about coupling diffusion, we now include the effect of deviations from ideality. Accounting for activity coefficients in the treatment has several consequences. It modifies the speciation of the metal cation under study, and therefore the diffusion rate of the various forms (free species and complexes). It also introduces a new driving force associated to the variation in space of the activity coefficients of the chemical species, and it modifies the magnitude of the internal diffusion electric field created by the major ions. The experimental observable is the result of these different effects.

In this work, we carried out experiments with the zinc(II) cation. This species has an important role in biological, geochemical and environmental media [START_REF] Pye | An ab initio Investigation of Zinc Chloro Complexes[END_REF][START_REF] Anderson | A microbeam XAFS study of aqueous chlorozinc complexing to 430 degrees C in fluid inclusions from the Knaumuehle granitic pegmatite, Saxonian granulite massif, Germany[END_REF]. It is well known for giving rise to multiple association with various anions [START_REF] Sillen | Stability Constants of Metal-Ion Complexes[END_REF] . Here, we used zinc chloride, which can give 4 complexes: ZnCl + , ZnCl 0 2 , ZnCl - 3 and ZnCl 2- 4 [START_REF] Sillen | Stability Constants of Metal-Ion Complexes[END_REF]. The zinc ion was taken in radioactive form, 65 Zn 2+ . A special technique consisting of a silica capillary inserted in a cylinder of scintillating plastic was utilized to observe the transient concentration profiles of the radioactive tracer [START_REF] Simonin | Closed capillary method for the diffusion of biological macromolecules[END_REF]. The major salt producing the internal electric field was LiCl. It was chosen because it produces a strong electric field, which is due to the greatly differing values of the diffusivities of the two ions Li + and Cl -. The observed experimental transient profiles were modeled by using finite difference (FD) and normal-mode (NM) methods. Deviations from ideality were included in the calculations by using the mean spherical approximation.

The main purpose of this paper is to investigate the validity of the modeling approach for the description of transient diffusion profiles for a multiply associating electrolyte. The theoretical descriptions include diffusion, ionic association and activity coefficients for all species. The maximum concentration of the major salt (LiCl) is sufficiently high to allow an appreciable influence of activity coefficients. However it is limited to 2M in order to minimize the effect of hydrodynamic interactions and other subtle effects such as reference frame or ionic relaxation [START_REF] Dufrêche | Analytical theories of transport in concentrated electrolyte solutions from the MSA[END_REF].

The remainder of this paper is organized as follows. The next section presents the experimental technique employed to probe the diffusion of the tracer in the cell, and the experimental conditions of this study. Then, a theoretical section exposes the basic equations governing the process and the methods used to calculate the observed effect. Next, the experimental data are compared with the results obtained from the model in the Results and Discussion section. Finally a conclusion summarizes the main results of this work and presents some prospects.

Experimental

Diffusion experiments on an ion which is a pure β emitter can be conducted adequately by using a cell made of perspex and scintillating plastic [START_REF] Simonin | Diffusion coupling in electrolyte solutions. 1. Transient effects on a tracer ion: sulfate[END_REF][START_REF] Simonin | Transport coupling of ions: Influence of ion pairing and pH gradient Application to the study of diagenetic fluxes[END_REF]. On the other hand if the tracer is both a β and γ emitter then it is necessary to stop the β's so that only the γ emission can be analyzed (otherwise it would be impossible to distinguish the respective effects of the β and γ radiations in the scintillating plastic). In the present work, the commercially available radioactive isotope of the ion studied here, 65 Zn 2+ , is precisely a β and γ emitter.

In order to solve this difficulty, experiments were carried out by using an adaptation of the closed capillary technique [START_REF] Simonin | Closed capillary method for the diffusion of biological macromolecules[END_REF]. It is depicted in Figure 1. The cell consists of a silica capillary (sealed at its bottom end) of 3 cm length and ca. 0.8 mm inner diameter which is introduced into a bore drilled through the center of a cylinder made of scintillating plastic. The plastic cell was designed in such a way that the points A, B and C are on the same straight line [START_REF] Simonin | Closed capillary method for the diffusion of biological macromolecules[END_REF]. The silica capillary (whose wall is ca. 1 mm thick) has the property of stopping the β + particles of 0.33 MeV energy emitted by 65 Zn 2+ . In this way, scintillation in the scintillating plastic is produced only by the γ radiation. In order to avoid capillarity problems at the mouth of the capillary, the latter was filled on a distance of 2.8 cm (i.e. 2 mm below the top of the capillary). Indeed, filling the capillary with liquid up to its mouth was observed to result in a significant loss of liquid in the course of the experiment. The top of the capillary was sealed with a piece of plastic paraffin film (Parafilm) to prevent evaporation. The silica capillary was adjusted within the plastic bore so that the mid-point of the liquid phase coincided with the top of the scintillating plastic cylinder.

The diffusion experiments were conducted as follows. The bottom part of the capillary (up to point O, see Figure 1), of length L, was filled with a solution containing 1M or 2M LiCl + 10 -2 M ZnCl 2 + radioactive zinc(II) in tracer amounts (of the order of 10 -8 M). Then, the top part of the capillary, of the same length L, was filled very carefully with a solution composed of 10 -2 M ZnCl 2 + radioactive zinc(II) at the same concentration. The solutions were introduced in the capillary by using a small plastic Pasteur pipette bearing a thin flexible Teflon tube at its end. This device allowed us to accurately handle very small volumes of liquid, of the order of 10 µL (half of the capillary volume). Next, the cell was introduced in a β radioactivity counter which determines the number of photons emitted in the scintillating plastic per Because it is expected that the γ rays can have only one efficient conversion in the plastic, it was shown [START_REF] Simonin | Closed capillary method for the diffusion of biological macromolecules[END_REF] that the probability of conversion of a γ ray is proportional to the length of plastic it crosses. This property allows an explicit calculation of the counting intensity as a function of the tracer concentration profile in the capillary. Then one can measure self-diffusion coefficients of γ-emitters with this technique (as was done before in the case of a radiolabelled biological molecule [START_REF] Simonin | Closed capillary method for the diffusion of biological macromolecules[END_REF]), and study diffusion processes like diffusion coupling.

The relative variation of the radioactivity intensity measured by the counter may be written as,

ρ(t) ≡ I(t) I(t = 0) -1 (1)
with I(t) the intensity measured by the photomultipliers at time t, so that ρ(t = 0) = 0. This intensity is produced by the radioactive tracer in its various forms. At a given position x and time t, we denote this total tracer concentration by C T (x, t).

With the notations of ref. [START_REF] Simonin | Closed capillary method for the diffusion of biological macromolecules[END_REF] we can write,

ρ(t) = ∫ 1 0 [θ -(z, t)f -(z) + θ + (z, t)f + (z)]dz ∫ 1 0 [f -(z) + f + (z)]dz (2)
in which z ≡ x/L for x > 0 (x-axis oriented upwards, see Figure 1) and z ≡ -x/L for x < 0, f -and f + are functions proportional to the 'efficacy' of a γ-ray emitted at position x (f -for x < 0 in the bottom part of the capillary and f + for x > 0 in the upper part), and θ -(z, t) and θ + (z, t) are the relative excesses of tracer concentration in bottom and upper parts of the capillary, respectively, that is

θ i (z, t) ≡ δC (i) T (z, t) C (0) T (3)
in which i = + or -and δC

(i) T ≡ C (i) T -C (0) T , C (-) T
and C

(+) T are the tracer concentrations for x < 0 and x > 0, respectively, and C (0) T is the initial uniform tracer concentration in the capillary. Expressions for the functions f -and f + were given in ref. [START_REF] Simonin | Closed capillary method for the diffusion of biological macromolecules[END_REF]. The concentration of the tracer is calculated below by using FD and NM techniques.

Theory

Basic relations in the non-ideal case

A mixture of two aqueous electrolyte solutions with common anion is considered. It is composed of a major salt, denoted by MA, constituted of monovalent M + and A -ions, and a salt TA n (T n+ cation with n ≥ 2) which is present in tracer amount in the solution. The major salt is totally dissociated. On the contrary, the T n+ cation is supposed to form N different complexes of the form TA

(n-k)+ k
, with k = 1, .., N , according to the following reactions,

T n+ + kA - TA (n-k)+ (β k ) (4)
characterized by the thermodynamic cumulative constants β k , with k going from 1 to N . Hereafter we will adopt the following notations for convenience. We will denote the free cation T n+ with subscript 0, the first complex TA (n-1)+ with subscript 1, ..., the k-th complex TA (n-k)+ k with subscript k. The subscripts + and -will be used for the cation and anion of the major salt, M + and A -, respectively.

Moreover, we will make the common assumption that the complexation reactions (eq. ( 4)) are very fast, so that a local equilibrium prevails for these reactions at any point in the capillary. Then, at any position x and time t, one can write for eq. ( 4),

C k = β k Γ k C 0 (C -) k (5)
in which any C i is the concentration of species i at x and t (C i (x, t)), and

Γ k ≡ γ 0 (γ -) k /γ k ( 6
)
where γ i stands for the activity coefficient of species i on molarity scale. In the experiments, the quantity of interest is the sum of the concentrations of the various forms of the tracer,

C T ≡ N ∑ k=0 C k (7)
In this work, besides deviations from ideality due to association itself, we take into account non-ideality effects on the transport of the species and on the formation of the complexes, arising from volume exclusion and electrostatic interactions,. This was done within the framework of the mean spherical approximation (MSA) as detailed below. Now we will write down the main equations governing the transport of the species in the solution.

The diffusion equations obeyed by the concentration C i of a species i is,

∂C i ∂t + ∇J i = Q i ( 8
)
with t the time, x the space coordinate shown in Figure 1, ∇ ≡ ∂/∂x the gradient operator in 1D, J i the flux density, and Q i the source term coming from the association reaction. In this equation, i denotes a free ion M + , A -, T n+ or a complex.

The flux density J i at position x and time t can be derived according to linear transport theory [START_REF] Turq | Normal Mode Analysis od Diffusion in Multicomponent Electrolyte Solutions. 1. Nonassociating Electrolytes[END_REF] to obtain the well known extended Nernst-Planck equation,

J i = -D i ∇C i -D i C i ∇ ln γ i + βD i C i z i eE (9)
in which D i is the diffusion coefficient of i, z i e is its charge, γ i is its activity coefficient, E is the local internal electric field and β = 1/kT with k the Boltzmann constant and T the temperature. The electric field is obtained from the dynamical electroneutrality condition for the major salt MA, J + = J - [START_REF] Miller | Application of Irreversible Thermodynamics to Electrolyte Solutions. I. Determination of Ionic Transport Coefficients l ij for Isothermal Vector Transport Processes in Binary Electrolyte Systems[END_REF] in which the fluxes of the minor species containing T are neglected (T in tracer amounts). At the same time the local electroneutrality condition is fulfilled, which gives

C + = C -= C s ( 11
)
again because one may neglect the concentration of the tracer, so that one gets by combining eqs. ( 9)-( 10),

E ≡ βeE = D + -D - D + + D - ∇ ln C s + 1 D + + D - ( D + ∂ ln γ + ∂C s -D - ∂ ln γ - ∂C s ) ∇C s (12
) for which we used eq. ( 11) and E has the dimension of the inverse of a length.

The combination of eqs. ( 9)-( 12) yields the well-known relation [START_REF] Harned | The Physical Chemistry of Electrolyte Solutions[END_REF][START_REF] Robinson | Electrolyte Solutions[END_REF],

J s = -D s ( 1 + ∂ ln γ s ∂ ln C s ) ∇C s ( 13
)
for the flux of the major salt s ≡ MA, with D s the Nernst-Hartley diffusion coefficient of the salt,

D s ≡ 2D + D - D + + D - ( 14
)
which is comprised between D -and D + . Eq. [START_REF] Turq | Normal Mode Analysis od Diffusion in Multicomponent Electrolyte Solutions. 1. Nonassociating Electrolytes[END_REF] shows that the effect of non-ideality is to modify locally the diffusion coefficient of the salt by a factor of (1+∂ ln γ s /∂ ln C s ).

The diffusion equation for the major salt is readily obtained from eq. ( 8) with Q + = Q -= 0 and eqs. ( 10), ( 11), [START_REF] Turq | Normal Mode Analysis od Diffusion in Multicomponent Electrolyte Solutions. 1. Nonassociating Electrolytes[END_REF],

∂C s ∂t = D s ∇ [( 1 + ∂ ln γ s ∂ ln C s ) ∇C s ] (15)
For the minor species containing the tracer we get using eqs. ( 8) and ( 9),

∂C i ∂t = D i ∆C i + D i ∇ [C i (∇ ln γ i -z i E)] + Q i ( 16
)
in which ∆ ≡ ∇.∇ = ∂ 2 /∂x 2 is the Laplace operator in 1D.

Non ideality was described using the mean spherical approximation (MSA) [START_REF] Simonin | Real Ionic in the Mean Spherical Approximation. 1. Simple Salts in the Primitive Model[END_REF], with values for the ion sizes and the solution permittivity that did not vary with salt concentration. A specific size was attributed to each complex (see the Results and Discussion section). The activity coefficients of the species were calculated at the McMillan-Mayer level (solvent of constant chemical potential and regarded as a continuum) and were converted to the experimental (Lewis-Randall) level by using the formula [START_REF] Simonin | Real Ionic in the Mean Spherical Approximation. 1. Simple Salts in the Primitive Model[END_REF][START_REF] Simonin | Study of experimental-to-McMillanMayer conversion of thermodynamic excess functions[END_REF][START_REF] Simonin | Real Ionic in the Mean Spherical Approximation. 2. Pure Strong Electrolytes up to Very High Concentrations, and Mixtures, in the Primitive Model[END_REF],

ln γ i = ln γ M M i -2C s Φ M M V i ( 17
)
with γ M M i and Φ M M the activity and osmotic coefficients, respectively, calculated at the McMillan-Mayer (MM) level within the MSA model, and V i the partial molar volume (PMV) of i. The PMVs of Zn 2+ and Cl -ions were estimated by taking the volume of the ion [START_REF] Huheey | Inorganic Chemistry: Principles of Structure and Reactivity[END_REF], which gives values of 0.001 and 0.012 L mol -1 , respectively. Indeed it was found that this simple choice gives the good order of magnitude for the PMV of ZnCl 2 at infinite dilution (0.0235 L mol -1 ). The PMVs of the complexes were computed from the formula,

V k = V Zn 2+ + kV Cl -.
Since the tracer is in very small amount in the solution, the activity coefficients of the complexes and of the salt ions only depend on the local concentration of the salt. Thus, in what follows, all physical quantities only depend on the local value of C s . Consequently the gradient of any quantity f may be calculated as, ∇f = f ′ × ∇C s [START_REF] Verweij | Free software developed by[END_REF] in which f ′ is the derivative of f w.r.t. C s , f ′ ≡ ∂f /∂C s . The diffusion equation for total tracer may be obtained from the equations for the various forms of the tracer. By considering that the concentrations of the complexes are equilibrated at any point in the capillary, we obtain by utilizing eqs. ( 5), ( 7) and [START_REF] Simonin | Diffusion coupling in electrolyte solutions. 1. Transient effects on a tracer ion: sulfate[END_REF],

C k = β k Γ k C s k C T / N ∑ n=0 β n Γ n C s n ( 19
)
In this relation we make the usual convention for n = 0, namely that β 0 = 1. We note that equilibration, expressed by eq. ( 19), is made very simple by the fact that T is in tracer amounts in the experiments. Summation of the diffusion equations for the free and associated forms of the tracer, together with the use of eqs. ( 12) and [START_REF] Sillen | Stability Constants of Metal-Ion Complexes[END_REF] in eq. ( 16), and also ∑ i Q i = 0, leads to the diffusion equation for total tracer. After some algebra one gets,

∂C T ∂t = D T ∆C T + d 1 ( ∇C s C s ) ∇C T + d 2 ( ∇C s C s ) 2 C T + d s ( ∆C s C s ) C T (20) in which D T = ∑ N k=0 D k β k Γ k C s k ∑ N k=0 β k Γ k C s k ( 21
)
is the (local) mean diffusion coefficient of the tracer, and the parameters d 1 , d 2 and d s have the dimension of a diffusion coefficient (but may be negative).

Their expression is given in the Appendix.

Normal-mode treatment

The normal-mode (NM) treatment consists of linearizing the transport equations. Any physical quantity, q(x, t), is decomposed in two terms: an equilibrium term, q ∞ (for t → ∞), and a non-equilibrium or perturbation term δq(x, t), that is, q(x, t) = q ∞ + δq(x, t). Linearization is implemented by keeping only the first-order terms in the equations.

Here, the concentrations of total tracer, C T (x, t), and of the salt, C s (x, t) are expressed as,

C i (x, t) = C ∞ i + δC i (x, t), in which C ∞ T = C (0) T and C ∞ s = C (0)
s /2 (with C (0) s the initial salt concentration in the bottom half of the capillary). The diffusion coefficients appearing in eq. ( 20) may be written as, D = D ∞ + δC s D ′ (with D = D T , d 1 , d 2 , d s ). Linearizing eq. ( 20) w.r.t. δC T and δC s leads to a diffusion equation for the total tracer in which the terms proportional to d 1 and d 2 disappear because they are of second order, δC s δC T or (δC s ) 2 , respectively (because ∇C i = ∇δC i for i = s, T ). So one obtains straightforwardly,

∂δC T ∂t = D ∞ T ∆δC T + d ∞ s C ∞ T C ∞ s ∆δC s (22)
in which D ∞ T and d ∞ s are the values of D T and d s (eqs. ( 21) and ( 32)) at equilibrium, for C s = C ∞ s . Analogously, the linearized diffusion equation for the salt is,

∂δC s ∂t = D ∞ s ∆δC s (23)
in which

D ∞ s ≡ D s ( 1 + ∂ ln γ s ∂ ln C s Cs=C ∞ s ) (24)
is the effective salt diffusivity at equilibrium. The coupled equations, eqs. ( 22) and ( 23), can be solved by using Fourier transform on space (denoted by an asterisk) and Laplace transform on time (denoted by a tilde), as exposed in previous work [START_REF] Simonin | Diffusion coupling in electrolyte solutions. 1. Transient effects on a tracer ion: sulfate[END_REF]. The result is,

δ C * T (q, s) = δ C * (s) T (q, s) + δ C * (m) T (q, s) (25)
with

δ C * (s) T (q, s) ≡ δC * T (q, t = 0) s + q 2 D * T (26) δ C * (m) T (q, s) ≡ -d ∞ s C ∞ T C ∞ s q 2 δC * s (q, t = 0) (s + q 2 D ∞ T )(s + q 2 D ∞ s ) (27)
in which q and s are the Fourier and Laplace variables associated to x and t, respectively, and δC * (s) T

is the self-diffusion term and δC * (m) T

is the migration term for the tracer.

We note that, in the present experiments, the initial tracer concentration profile is uniform (equal concentration in the two halves of the capillary), so that δC T (x, t = 0) = 0 and thus δC * T (q, t = 0) = 0. Therefore, only the migration term, δ C * (m) T , will remain in Eq. ( 25). Finally, by using the same method as in previous work [START_REF] Simonin | Diffusion coupling in electrolyte solutions. 1. Transient effects on a tracer ion: sulfate[END_REF], we find that the solution for the function ρ(t) (Eq. ( 1) or ( 2)) is,

ρ N M (t) = - 4 π d ∞ s ∞ ∑ k=0 F k 2k + 1 exp(-q 2 2k+1 D ∞ T t) -exp(-q 2 2k+1 D ∞ s t) D ∞ s -D ∞ T , (28)
where q n ≡ nπ/2L, and

F k = ∫ 1 0 [f -(z) -f + (z)] sin [ (2k + 1) π 2 z ] dz / ∫ 1 0 [f -(z) + f + (z)] dz (29)
In eq. (28) d ∞ s was computed from eq. ( 32) at equilibrium and the quantity Ψ ′ k (appearing in eq. 32 for d s ) was calculated by numerical differentiation.

For the geometry of the set-up depicted in Figure 1, the coefficients F k were found to have the following values [START_REF] Simonin | Closed capillary method for the diffusion of biological macromolecules[END_REF] : F 0 = 0.472677, F 1 = 0.0564312, F 2 = 0.0233391, F 3 = 0.00424871, F 4 = 0.004660106, F 5 = 0.0000347272. Since all these parameters are positive, we deduce from eq. ( 28) that ρ N M has the sign of -d ∞ s . For sufficiently long times, eq. ( 28) may be reduced to its first exponential term. In the present case we find that this approximation is valid to 1% accuracy for t 10 h. By differentiation of this simplified expression w.r.t. time one gets the order of magnitude of the time t m at which ρ N M is maximum (in absolute value),

t m = 4L 2 π 2 D ∞ s ln x x -1 (30) with x ≡ D ∞ T /D ∞ s .
By inserting this expression for t m into eq. ( 28) we find after some simplifications that the maximum amplitude of ρ N M (at t = t m ) is approximately given by,

ρ N M (t = t m ) ≃ - 4 π F 0 d ∞ s D ∞ s x x/(1-x) ( 31
)
The function x x/(1-x) decreases monotonously with x from a value of 1 for x = 0 to 0 for large values of x.

Finite difference method

A classic finite difference algorithm was used in this work. The main ingredients of the method were described in our previous study [START_REF] Simonin | Diffusion coupling in electrolyte solutions. 1. Transient effects on a tracer ion: sulfate[END_REF]. The diffusion equations were discretized in space and time, and the system was symmetrized at both ends of the capillary. A number of N c = 200 cells was taken to represent the capillary. The result was observed to be nearly independent of N c for N c ≥ 200. At each time step, all species were let to diffuse individually according to eqs. ( 15) and ( 16) with Q i = 0, and then the concentrations of the complexes were equilibrated in each cell.

A difference in the present work comes from the inclusion of non-ideality. In our earlier treatment [START_REF] Simonin | Diffusion coupling in electrolyte solutions. 1. Transient effects on a tracer ion: sulfate[END_REF], conservation of matter in the capillary was ensured by imposing a zero electric field at both ends of the capillary. In the present work, it was found that this condition could be satisfied by discretizing eq. ( 16) (without expanding it) and by imposing the relation, ∇ ln γ i -z i E = 0, in cells no. 1 and N c at every time step.

Results and Discussion

Parameter values

First, parameter values had to be determined for the description of association, activity coefficients and transport of species in the model. An approximate procedure had to be used in order to reduce the number of parameters in the model. Their values are collected in Table 1. They were determined as follows.

The diameters of the salt ions Li + and Cl -were required for the computation of activity coefficients of all species within the MSA [START_REF] Simonin | Real Ionic in the Mean Spherical Approximation. 2. Pure Strong Electrolytes up to Very High Concentrations, and Mixtures, in the Primitive Model[END_REF]. In this framework, we usually keep the diameter of Cl -ion equal to its Pauling crystallographic value (3.62 Å) because it is essentially unhydrated. Then, the diameter of (hydrated) Li + ion was adjusted by fitting the experimental osmotic coefficient of aqueous LiCl solution up to 2M, by keeping the dielectric constant of solution equal to that of water.

The values for the cumulative constants β n were taken from a database published in the geochemical community for the determination of speciation in natural waters [START_REF] Turner | The equilibrium speciation of dissolved components in freshwater and seawater at 25 • C and 1 atm pressure[END_REF]. These values are close to those proposed by Marcus [START_REF] Marcus | Anion Exchange of Metal Complexes. VIII. The Effect of the Secondary Cation. The Zinc-Chloride System[END_REF].

A fit of experimental osmotic coefficient data for zinc perchlorate solution up to 0.2 mol kg -1 [START_REF] Goldberg | Evaluated Activity and Osmotic Coefficients for Aqueous Solutions: Bi-Univalent Compounds of Zinc, Cadmium, and Ethylene Bis(Trimethylammonium) Chloride and lodide[END_REF] within the MSA with constant solution permittivity gave an optimum zinc(II) ion diameter of ∼ 6.4 Å with a diameter of 4.53 Å for the perchlorate ion [START_REF] Simonin | Real Ionic in the Mean Spherical Approximation. 2. Pure Strong Electrolytes up to Very High Concentrations, and Mixtures, in the Primitive Model[END_REF]. In this fit, zinc perchlorate was regarded as being not associated as found experimentally [START_REF] Jones | A Search for Perchlorate Complexes. Raman Spectra of Perchlorate Solutions[END_REF][START_REF] Ohtaki | X-Ray Diffraction Studies of the Structure of Hydrated Divalent Transition-Metal Ions in Aqueous Solution[END_REF]. This value of 6.4 Å was taken here for the computation of the activity coefficient of free zinc(II) ion in LiCl solution.

It is known from experiment that the structure of the zinc(II) ion in water is a hexa-coordinated aquo-complex, ZnW 2+ 6 [START_REF] Ohtaki | X-Ray Diffraction Studies of the Structure of Hydrated Divalent Transition-Metal Ions in Aqueous Solution[END_REF] (in which W represents the water molecule), with an average lifetime of a few nanoseconds [START_REF] Salmon | The dynamics of aqueous Zn 2+ solutions: a study using incoherent quasi-elastic neutron scattering[END_REF]. The maximum diameter of this species is of the order of 1.76 + 2 × 2.76 ∼ 7.3 Å, as found by adding the diameters of the bare zinc(II) ion [START_REF] Huheey | Inorganic Chemistry: Principles of Structure and Reactivity[END_REF] and two times that of the water molecule. Besides, the diameter of a sphere having the same volume as this species is ∼ 5.1 Å. We note that the MSA adjusted diameter of the zinc(II) ion in perchlorate solution is nearly the average of these two values. So, the diameter of 6.4 Å may be interpreted as the effective size of the aquo-complex in the MSA.

The structure of zinc chloride complexes in water has been studied experimentally [START_REF] Maeda | The Structure of Zinc Chloride Complexes in Aqueous Solution[END_REF] and by ab initio numerical simulation [START_REF] Parchment | Speciation in Aqueous Zinc Chloride. An ab Initio Hybrid Microsolvation/Continuum Approach[END_REF][START_REF] Pye | An ab initio Investigation of Zinc Chloro Complexes[END_REF][START_REF] Harris | Zinc Complexation in Hydrothermal Chloride Brines: Results from ab Initio Molecular Dynamics Calculations[END_REF]. The consis-tent picture emerging from these studies is one in which the monochloro and the dichloro complexes are hexa-coordinated, and the trichloro and tetrachloro complexes are tetra-coordinated. Their formulas are ZnClW + 5 , ZnCl 2 W 0 4 , ZnCl 3 W -and ZnCl 2- 4 . The change of coordination in the complexes with 3 and 4 chlorides is likely to originate from the fact that it lowers the electrostatic repulsion between the charges on the chlorides as compared to hexa-coordinated compounds [START_REF] Harris | Zinc Complexation in Hydrothermal Chloride Brines: Results from ab Initio Molecular Dynamics Calculations[END_REF]. The diameters taken for the complexes in the MSA model were estimated from the MSA diameter of zinc(II) ion (σ 0 Å) and the volumes of the complexes (v k ) by using the formula,

σ k = σ 0 (v k /v 0 ) 1/3
. This led to the values listed in Table 1.

The diffusivities of Li + , Cl -and Zn 2+ at 25 • C were obtained from the values for the limiting equivalent conductivities [START_REF] Robinson | Electrolyte Solutions[END_REF] by using a classic conversion formula (eq. 11.49 of ref. [START_REF] Robinson | Electrolyte Solutions[END_REF]). Those for the complexes were assessed by utilizing the D value for the zinc(II) ion, the values of the complex diameters (σ k ), and the Stokes-Einstein formula, from which we got, D k = D 0 σ 0 /σ k , by assuming that σ k is a good approximation for the hydrodynamic diameter of the complex k. It is noticed in Table 1 that the value of D is not monotonous in the series because of the transition from hexa-to tetra-coordination when passing from ZnCl 2 to ZnCl - 3 .

Forces in play

In order to gain more insight into the process, the speciation found for the zinc ion at LiCl concentrations of 0.5M, 1M and 2M from the MSA description, and in the ideal case with association is shown in Table 2. The first two concentrations, 0.5M and 1M, correspond to the values for the concentration of LiCl at point O (see Figure 1) when the bottom of the capillary is filled with 1M or 2M LiCl solution at t = 0. Indeed the LiCl concentration at point O is approximately constant in time during the diffusion process and it is equal to half of the initial concentration in the bottom part of the capillary. It must be noticed that this zone of the capillary is where most of the internal electric field effect is created. This is due to the fact that the electric field is the most intense in this region of highest gradient of LiCl. This field originates from the diffusion of the major ions. The chloride ion diffusivity being larger than that of the lithium ion, an electric field is created, directed from the bottom to the top of the capillary, that pulls the slow ion (the positive lithium ion) and slows down the fast ion (the chloride ion).

It is seen in this table that as expected the speciation of zinc is shifted towards the formation of complexes of higher stoichiometry when C LiCl is increased. In the upper region of the capillary, as long as LiCl has not reached this zone, the zinc is in the form of Zn 2+ ion (dilute ZnCl 2 solution). It is seen also in this table that the effect of deviations from ideality is rather weak at 0.5M LiCl, moderate at 1M and appreciable at 2M.

The value of the average valency of zinc containing species, ⟨z⟩, is given in the last column of Table 2. It has a value of 2 at very low C LiCl , at which the free zinc(II) ion is the dominant form of zinc. This quantity gives an indication of the net effect of the electric field on the zinc species (taken globally, in all forms) because the complexes are expected to have similar diffusivities. In particular, since ⟨z⟩ is positive at 0.5M LiCl and negative at 1M, we find that the electric field should have opposite effects on zinc at these two concentrations (which correspond to those at point O in the capillary for initial 1M and 2M LiCl, respectively, in the bottom part). At very low LiCl concentration, the zinc(II) ion would be dragged in the direction of the electric field (from the bottom to the top of the capillary) because it is positive, which would make the function ρ take negative values. When C LiCl is increased, ⟨z⟩ is also decreased. Then we expect that ρ will take less negative values for initial 2M LiCl than for 1M under the effect of the electric field.

Besides the influence of the electric field, the second main phenomenon governing the motion of zinc is self-diffusion (the term -D i ∇C i in eq. ( 9)). In solutions of 1M and 2M LiCl initially in the bottom of the capillary, the Zn 2+ ion is the minor zinc species and it is by far the major ion in the top part. This strong gradient causes natural diffusion of Zn 2+ ion from the top to the bottom half of the capillary. It is the reverse situation as regards the complexes that are nearly inexistent in the top part. This induces crossed fluxes of free zinc(II) ion and of the complexes. However the overall flux of zinc from the bottom to the top of the capillary is likely to be small because these fluxes may compensate each other.

The last effect is that of deviations from ideality on the fluxes as expressed in eq. ( 9) by the term -D i C i ∇ ln γ i . The activity coefficients of the zinc compounds obtained from the MSA model are collected in Table 3. The bivalent ions Zn 2+ and ZnCl 2- 4 have the lower values because they have higher electrostatic interactions with the major salt, which gives them lower electrostatic energies. The monovalent ions ZnCl + and ZnCl - 3 , and still more the neutral ion ZnCl 2 , have higher activity coefficient values. According to eq. ( 9) the flux of species i due to non-ideality has the sign of -∇ ln γ i . Thus, with the x-axis being oriented from the bottom to the top of the capillary (see Figure 1), the sign of the overall flux for a species i may be expected to have basically the sign of ln γ i reported in Table 3. Consequently, according to this table, we see that the species should behave differently with respect to this effect. On the other hand, we observe that the values of ln γ i all increase in algebraic value from 1M to 2M LiCl concentration. Therefore the effect of activity coefficients should be stronger at 2M LiCl than at 1M. The transfer of zinc under this effect from the bottom to the top of the capillary should be larger at 2M than at 1M , so making ρ more negative.

This preliminary discussion shows that the transport of zinc in the capillary is the result of several antagonistic effects. At low C LiCl , the main effect would be due to the electric field (which remains finite by virtue of eq. ( 12) and becomes independent of C s when C s becomes very small), so that ρ would be negative. When C LiCl is increased, the effect of self-diffusion must be small, that of the electric field is to reduce the amplitude of ρ, and that of activity coefficients is to enhance it. The theoretical FD and NM treatments were utilized to solve this problem. The results are exposed in the next section.

Comparison of experimental and theoretical results

Typical experimental curves (denoted by the letter "e") and theoretical results are presented in Figures 2 and3, in which -ρ is plotted as a function of time. The fact that the experimental plots appear as broken lines originates from the natural fluctuations in γ-emission by the radioactive zinc tracer. In Figure 2, the FD result including association and deviations from ideality is in very good agreement with the experimental data. The FD curves obtained in the ideal case (FD ideal), and without association (FD ass0), are also shown. They indicate how sensitive the FD result is to the incorporation of these two phenomena. It is noticed that the ideal plot for -ρ is below the non-ideal result (maxima of ca. 6% and 8%, respectively). On the contrary, in the absence of association, the result (FD ass0) is significantly above the FD plot, mainly because associated forms of zinc ion have a lower (possibly negative) electric charge than the free zinc(II) ion. Consequently the effect of the electric field is lower on these species or it is in the reverse direction. The normal modes (NM) plot is slightly above the FD plot, with a relative difference of ∼ 17% between the two maxima. Moreover, the positions in time of the maxima for -ρ are nearly identical (∼ 22.0 h and 22.6 h for FD and NM, respectively), and the value for NM coincides with the one given by eq. ( 30) (t m = 22.6 h). The value for the NM maximum is accurately given by eq. [START_REF] Huheey | Inorganic Chemistry: Principles of Structure and Reactivity[END_REF].

In Figure 3, the same curves as in Figure 2 are plotted in the case of initial 2M LiCl. It is observed that the full FD result (accounting for association and activity coefficients) is close to the experimental curve ("e"). At the same time, the sensitivity of the FD result to the inclusion of association and deviations from ideality is much more pronounced than at 1M LiCl (see Figure 2). Here, the FD result in the ideal case is much below the FD plot, which shows the importance of including non-ideality in the description. Moreover, the formation of zinc complexes is seen to cause a large drop in the observed maximum for -ρ. In the absence of association we find that the maximum for the FD plot is of ∼ 18% as compared to ∼ 10% for the experiment and the full FD curve. Thus, the full FD calculation including association and activity coefficients seems capable of capturing the main effects acting on the various forms of the solute. Besides, the NM plot is just a bit higher than the FD plot, with a relative difference between the two maxima of ∼14%. The positions in time of the maxima are a bit more different from each other than at 1M LiCl (∼ 17 h and 19 h for FD and NM, respectively). This value for NM is precisely given by eq. ( 30) (t m = 19.2 h).

Conclusion

It has been found that the model developed in this work provides a good prediction of transient diffusion profiles for a multiply associating ion forming chemical complexes in a supporting electrolyte. This satisfactory description was developed by taking literature values for the association constants and by describing deviations from ideality within the MSA with reasonable input parameter values. The prediction consisted of using a finite difference algorithm and, besides, an approximate normal mode analytic calculation. Both accounted for diffusion, internal electric field and activity coefficients for all species. The normal mode result was found to be a little larger than the more accurate finite difference result. The NM method may be used to estimate the order of magnitude of the experimental effect and it provides a confirmation of the FD result.

In subsequent work, we will study the case of electrolytes forming several ionic complexes. In such systems, the cation interacts strongly with an anion through electrostatic forces but does not form covalent complexes. This is for instance the case of lanthanides with some halides. It will be interesting to examine whether coupling diffusion can be described within the classic chemical model [START_REF] Barthel | Physical Chemistry of Electrolyte Solutions, Modern Aspects[END_REF] in which an ionic complex is regarded as a true chemical complex.

Figure 1 :

 1 Figure 1: Sketch of the diffusion cell. The points A, B and C are aligned.

Figure 2 :

 2 Figure 2: Plot of -ρ as a function of time for the experiment with 1M LiCl. (e) Experimental points; (FD) and (NM) FD and NM results with association and deviations from ideality, respectively; (FD ass0) FD without association but with deviations from ideality; (FD ideal) FD for ideal solution with association.

Figure 3 :

 3 Figure 3: Plot of -ρ as a function of time for the experiment with 2M LiCl. (e) Experimental points; (FD) and (NM) FD and NM results with association and deviations from ideality, respectively; (FD ass0) FD without association but with deviations from ideality; (FD ideal) FD for ideal solution with association.

Table 1 :

 1 Values taken for diameters (σ i ), diffusivities (D i ), and cumulative formation constants β i of species.

	Li +				4.26	1.03 a
	Cl -				3.62	2.03 a
	ZnW 6	2+			6.40	0.703 a
	ZnClW 5	+	6.65	0.677	0.49
	ZnCl 2 W 4	0	6.88	0.654	0.62
	ZnCl 3 W -	6.50	0.692	0.51
	ZnCl 2-4			6.75	0.667	0.20

Species

σ i ( Å) D i (10 -5 cm 2 s -1 ) log 10 β i a Experimental value

Table 2 :

 2 Speciation and average valency for zinc at various concentrations of LiCl, obtained from the present MSA model (first line for every LiCl concentration) and in the ideal case (second line, in parentheses). Average valency of zinc species in its various forms

	C LiCl	Zn 2+	ZnCl +	ZnCl 2	ZnCl -3	ZnCl 2-4	⟨z⟩ a
	0.5M 25.4%	40.8%	19.3%	11.6%	2.9%	0.74
		(23.2%) (37.3%) (26.2%) (10.6%) (2.7%) (0.68)
	1M	10.5%	24.0%	19.8%	27.7%	18.0%	-0.19
		(7.3%) (23.1%) (31.8%) (25.2%) (12.6%) (-0.13)
	2M	2.4%	4.4%	6.5%	22.2%	64.5%	-1.4
		(1.3%) (8.1%)	(22%) (34.5%) (34.1%) (-0.92)

a

Table 3 :

 3 Values of ln γ i for the zinc compounds obtained from the present model.[LiCl] Zn 2+ ZnCl + ZnCl 2 ZnCl -

	3	ZnCl 2-4

in the case of Mg 2+ and SO2- 4 [11], or to the hydrogen sulfate ion HSO - 4 in the case of sulfate in the pH gradient.

Ackowledgement

This work was supported by a grant from ECOS-Nord (project M11P03).

Appendix

The coefficients appearing in eq. ( 20) are as follows,

where a prime denotes differentiation w.r.t. C s , i.e. f ′ ≡ ∂f /∂C s , and we have introduced the following definitions

Note that, in these relations,

The derivatives in these equations involve differentiation of individual activity coefficients w.r.t. the salt concentration C s . They may be calculated for a given model of deviations from ideality, either analytically (if possible) or numerically.


Wael Ben Messaoud
email: wael.ben-messaoud@ed.univ-lille1.fr

Michel Amberg

Betty Lemaire-Semail

Frédéric Giraud

Marie-Ange Bueno

High fidelity closed loop controlled friction in SMARTTAC tactile stimulator

Keywords: «Tactile feedback», «Force Control», «Closed-Loop Control», «Sliding mode control», «Piezo actuators»

The present study investigates the implementation of a closed loop control in a new smart tactile stimulator called "SMARTTAC". This haptic device is able to stimulate the fingertip in order to simulate the texture perception of different real surfaces such as gratings. Two closed loop control are implemented to ensure both a stable vibration amplitude despite of external disturbances and accurate control of the friction coefficient using the new design of the tactile stimulator. Experimental measurements give satisfactory results both on dynamics and robustness points of view. High fidelity closed loop controlled friction in SMARTTAC tactile stimulator BEN MESSAOUD Wael EPE'15 ECCE Europe

Introduction

Nowadays, the sense of touch has been increasingly explored by the design of tactile devices aiming to render the sensation of touch of real surfaces using a programmable device. Different characteristics of the surface can be reproduced in the experience of touch, for example temperature, softness or texture friction characteristics. Three major types of friction modulation devices can be classified depending on the operating technology to produce the tactile sensation. The surface acoustic wave is a technology able to change the fineness of the surface's grain [START_REF] Nara | Surface Acoustic Wave (SAW) tactile display based on properties of mechanoreceptors[END_REF], the electrovibration technology increases the friction by creating an attractive electrostatic force between the sliding finger and the surface [START_REF] Bau | TeslaTouch: electrovibration for touch surfaces[END_REF], [START_REF] Vezzoli | Electrovibration modelling analysis[END_REF] and the ultrasonic friction-modulation devices operates to reduce the friction by creating an air film lubricating the contact finger and the surface [START_REF] Winter | Optimal design of a squeeze film actuator for friction feedback[END_REF], [START_REF] Watanabe | A method for controlling tactile sensation of surface roughness using ultrasonic vibration[END_REF]. The new tactile stimulator called "SMARTTAC" is a device able to stimulate the fingertip exploring its surface in order to simulate different textures. This device is an ultrasonic friction modulation device because it is based on producing ultrasonic vibration of a plate above which the finger is sliding. The objective is to use this device for fabrics simulation, particularly for the e-commerce shopping. One of the significant difference between reality and virtual shopping is that the purchaser cannot have the opportunity to touch and compare products [START_REF] Lee | Behavioral Model of Online Purchasers in E-Commerce Environment[END_REF]. In that case, a virtual tactile feedback provided by a tactile stimulator could be helpful. But to be able to simulate such complex textures, the tactile feedback device must be accurate and reliable. So the aim of this study is to propose a control strategy in order to cope with the numerous phenomena which may disturb the touch feedback. Such a stimulator could be used also to help people who suffer from tactile disorders for diagnosis and rehabilitation processes. On a technical point of view, this tactile stimulator is based on the creation of an air gap between the finger a vibration of the surface. This vib fulfil the requirement of high fi assuming the vibration amplitude to cope with external disturbance static performances, and reducti measure the friction and norm comparison to an external tribom finger-to-stimulator friction. In this paper, the new tactile stim will be exposed before proceed controlled system will be carrie control and on the other side, the

Description of the "SMA

The SMARTTAC is a develop piezoelectric ceramics to make piezoceramics (four in each extr stimulator is sized 76×41×1.2 m smart tactile stimulator because lateral force between the finger lateral force measurement, two placed in a compact and deform placed under the plate (Fig. 1b). calibration the tangential force ap by a support frame (A), tighten elastically connected to a central This elasticity is allowed thanks sliding. The normal force applied vibration nodes (S3-S6). The fri law on each sampling period. To stimulator based on the friction r (a) Fig. 1 a): Photo of the SMARTTA vibrating plate, able to measure u it explores the active surface of " nd a flat vibrating plate, thanks to a high frequ bration is able to reduce the friction coefficient [4 fidelity simulation, the device is equipped with e control as well as the tuning of the operating fr es during the exploration process, while keeping a ion of energy consumption as well [START_REF] Ben Messaoud | Closed-Loop Control for Squeeze Film Effect in Tactile Stimulator[END_REF]. Moreove mal forces during tactile task, as a tribometer meter is that it is possible to control the vibration mulator will be presented, then, the modelling an ding to the control of the system. Finally, expe ed out showing on one side the accuracy of the e first attempts for actual friction regulation.

ARTTAC" stimulator

ped version of the STIMTAC stimulator [10 a plate vibrating at one of its ultrasonic reson remity) are glued under the active surface of the mm 3 (Fig. 1a). In comparison with STIMTAC, it integrates new sensors allowing a real time rtip and the vibrating plate during the explorat strain gauge sensors (S1) and (S2) (FSS1500 f mable mechanical structure, especially designed S1 and S2 measure respectively the forces in x an pplied on the plate is obtained. The mechanical s ned by two screws to the top and bottom of t l part (C) attached by two more screws to a solid s to six deformable slats (B) between (A) and (C d by the finger is measured from four other strain iction coefficient (COF) can be easily deduced o sum up, the "SMARTTAC" is the unique devi eduction and a tribometer.

(b)

AC, b): Fabricated mechanical structure placed u using two force sensors, the lateral force applied b "SMARTTAC" uency, low amplitude 4], [START_REF] Biet | Implementation of tactile feedback by modifying the perceived friction[END_REF], [START_REF] Winfield | T-PaD: Tactile Pattern Display through Variable Friction Reduction[END_REF]. In order to closed control loops frequency; this allows accurate dynamic and er, SMARTTAC can r. The advantage in n, as a function of the nd identifying process erimental tries of the e vibration amplitude 0]; it is actuated by nant frequency. Eight e aluminum plate; the "SMARTTAC" is a measurement of the tion process. For the from Honeywell) are d for this application, nd y directions. After structure is composed the tactile plate, and d aluminum base (D). C) when the finger is gauges located at the thanks to Coulomb's ce which combines a under the active by the fingertip when

Vibration amplitude modelling

In this section, we propose a modelling of the vibration amplitude behaviour to compensate the disturbance factors influencing on the vibration of the stimulator such as the influence of the finger pressure and the influence of the resonant frequency due to the variation of temperature [START_REF] Ben Messaoud | Closed-Loop Control for Squeeze Film Effect in Tactile Stimulator[END_REF], [START_REF] Wiertlewski | A High-Fidelity Surface-Haptic Device for Texture Rendering on Bare Finger[END_REF]. So, the system can be seen as Multiple-Input Multiple-Output or MIMO. The two inputs are the components of the voltage supply, which are the voltage amplitude and the instantaneous excitation angular frequency . In other words, the system can be represented in a rotating d-q frame, with two components and ; this representation is classical for the study of electromagnetic machines [START_REF] Leonhard | Control of Electrical Drives[END_REF] and it was also operated in the case of a standing wave piezoelectric actuator in [START_REF] Pigache | Modelling and identification of a planar standing wave ultrasonic motor[END_REF]. The two outputs are the two components of the vibration amplitude which are the real and the imaginary parts respectively and . The model of the tactile plate associated with a vibration mode is given by the following second order equation [START_REF] Pigache | Modelling and identification of a planar standing wave ultrasonic motor[END_REF]: [START_REF] Nara | Surface Acoustic Wave (SAW) tactile display based on properties of mechanoreceptors[END_REF] where is the modal mass, is the modal damping, is the electromechanical conversion factor and represents the modal stiffness. The vibration amplitude and the voltage can be represented with their complex form as following (the d-q axes are not yet defined):

(

To simplify the modelling and the control, we choose the d-q frame attached to the voltage supply vector, such as 0. By replacing the complex format of and ( 2) and ( 3) in (1), the real and imaginary parts of the equation ( 1) can be written as follows:

2 0 (4) 2 (5)
At the resonance, 0 because a 90° is between the vector and [START_REF] Giraud | Design and control of a haptic knob[END_REF].The modelling process of the system of equation [START_REF] Watanabe | A method for controlling tactile sensation of surface roughness using ultrasonic vibration[END_REF] gives the transfer function between and which is a first order transfer function represented as follows:

⁄ 1 2 ⁄ (6)
with the angular resonant frequency. At vicinity of the resonance, if we control the based on the previous equation ( 6) so as to maintain it constant, their derivative terms of are zero. So, we extract from the equation (4) the model between / as a function of the variation of the angular frequency by considering , synthetized as follows:

2 ⁄ 1 2 ⁄ (7)
From the two models, we observe that the constant time of the two first order equations is the same. It may be noted also that the gain of the transfer function ( 7) is equal to the negative value of the constant time. It is the major advantage of this frame choice to give such a simple modelling, easy to identify experimentally. After the identification of the process modelling, a vector control of has been performed, leading to a high dynamic performance system and an accurate robustness to face finger pressure and temperature variation disturbances. Experimental results are given later to assess to performances of the controllers. Nevertheless, these performances do not guarantee the accurate stimulation of the user, which highly depends on the finger-plate interaction and may vary from one user to the other. Another additional control must then be performed.

Friction coefficient control Problematic

After tribological measurement of the coefficient of friction (COF) with 5 subjects (aged from 26 to 28, average of 27) at 5 controlled vibration amplitudes [0; 0.5; 1; 1.5; 2] µm, it can be observed in Fig. 2 that the percentage of the ratio between the reduced friction μ and the initial COF μ is highly variable depending on each person and on the conditions of measurement (the moisture for example increases the friction [START_REF] Pasumarty | Friction of the Human Finger Pad: Influence of Moisture, Occlusion and Velocity[END_REF], [START_REF] Tomlinson | Understanding the Friction Mechanisms Between the Human Finger and Flat Contacting Surfaces in Moist Conditions[END_REF]). The initial friction represents the friction between the finger of the user and the surface when the vibration is switched off whereas the reduced friction is the friction when the vibration is switched on at different level from 0 to 2 µm with a step of 0.5 µm.

Fig. 2:

The ratio between the reduced friction μ and the initial friction μ in % relative to the vibration amplitude for five subjects plotted by colored points. The mean of measurements of all subjects is plotted by a continuous red curve.

Our aim is to propose a tactile stimulator able to render accurately texture feelings whatever the touch conditions or the user finger characteristics. So, we propose a direct control of the COF in order to standardise the friction reduction which may be provided by the tactile feedback device. Using this smart stimulator called SMARTTAC described in the first part of this paper; it is possible to measure the friction and the normal forces leading to a real time measurement of the instantaneous friction coefficient μ / between the finger and the surface of the plate. By knowing the instantaneous friction under the finger, an adjustment of the vibration amplitude can be achieved to impose a desired friction under the finger.

Control method

Because of the non-linear gain of the system having the vibration amplitude as input and the friction as output, it was not possible to implement a standard controller. So, we develop a controller based on data measurements. The sliding mode control is operated due its capability to operate with nonlinear systems [START_REF] Shtessel | Introduction: Intuitive Theory of Sliding Mode Control[END_REF]. It is placed in a closed loop including the previous vibration amplitude control (see Fig. 4). The controller acts on the vibration amplitude to control the friction coefficient between the finger and the plate. The reference of the controller is the difference between the instantaneous and the desired frictions divided by the desired friction coefficient as:

100 μ μ μ (8)
The controller task is to minimise the parameter to impose the instantaneous desired friction. The controlled variable by the controller is the parameter ∆ , it represents the value to add to the vibration amplitude to reach the friction reference. The developed controller operates at three cases; when 25%, the controller's output is 1 µm, if this ratio is up to 25%, the controller output is 1 µm. the last case is when the ratio is between -25% and 25%, in this case the controller evolves according to the equation 0.04 . The controller output is limited by a saturation set at ±1 µm to avoid a possible peak of vibration under the finger. The constant 0.04 multiplied by the parameter was chosen is order to find a compromise between the minimisation of the response time and the reduction of the chattering problem when the plate was controlled [START_REF] Utkin | Chattering Problem in Sliding Mode Control Systems[END_REF]. To sum up, the controller output is the following:

∆ 1, 25 0.04 , 25 25 1, 25 (9)
The figure 3 illustrates the controller output which is the vibration amplitude as a function of E which is the controller reference and can be considered as the friction error.

Fig. 3: Controller output acting on the vibration amplitude depending on the friction error

To sum up, the closed loop control strategy method is based on two nested loops. The first one consists on controlling the vibration amplitude using the controller which acts on the voltage applied to actuators. The control method is based on modelling the system in a d-q frame to determine the transfer functions between inputs and outputs. The goal is to maintain constant and equal to the desired vibration amplitude whatever the disturbances and to vary the excitation frequency to guarantee the system operation at resonance. The main loop is based on a direct control of the friction between the sliding finger and the plate. The controller acts on the desired vibration to ensure a desired friction coefficient μ μ.

Control implementation

This section aims at giving expe second on the friction coefficient

Vibration amplitude control

Based on the identification of th controller adjusts the voltage in second controller is able to cont resonance which maximizes the together at 1 kHz frequency realperformances and robustness. As control law against the temperatu to temperature variations which frequency is a real problem in us excited at the resonance and th variations have been described in voltage at an initial frequency although the real temperature dur frequency). The closed loop sh minimum time and as a conseq excitation voltage. Fig. 5 illustr amplitude and angular frequency a second experiment for which whereas we work at 20°C), and amplitude for the two experiment ategy aiming to a control of the vibration amplitud ariation of the resonance frequency and the influe splacement when the finger is sliding on the plate, n amplitude to control the friction.

and experimental results

erimental results obtained first for the vibration a t control. l he two transfer functions, two PI controllers were n order to obtain a constant vibration amplitude rol the parameter to guarantee the excitation e ratio vibration/voltage (q axis). The controll -time. The PI controller parameters are chosen to s an example here, we have chosen to illustrate ure variation: as a mechanical system, the tactile s induce a shift in the resonant frequency. A vari se if it is not taken into account. In fact, the active he reduction of the friction is not at its maxim n [START_REF] Vezzoli | Electrovibration modelling analysis[END_REF]. To check the robustness of our control law, which corresponds to the resonant frequency ring this experiment is about 20°C (that means a hould lead to a tracking of the actual resonan quence, the desired vibration amplitude will b rates the fully controlled system behaviour wi y (red curve). The second result in Fig. 5 (green c h the system is excited at the resonant frequen d only the vibration amplitude is controlled. T ts was 1 µm PtP. de in the internal nce of the finger , while the main amplitude control and e calculated. The first e (d axis), while the n of the system at the ers are implemented ensure both dynamic the robustness of the stimulator is sensitive iation of the resonant e surface is no longer mum potential. These , the idea is to supply at 30°C (43250 Hz), 0.3% higher resonant nce (43390 Hz) in a e obtained with less th the two loops, in curve) corresponds to ncy 43250 Hz (30°C The desired vibration Fig. 5: Experimental comparison between vibration amplitude response of the controlled vibration system, and the response in the case of the double controlled system against temperature variation and vibration amplitude damping.

In static, it can be observed that in both experiments, the vibration amplitude is maintained constant (1 µm). But, in dynamic, we can observe a difference in the 5% response time of the system (5 ms in the vibration controlled system against 2 ms for the system controlled in vibration and frequency). To evaluate the importance of the double control in the energy consumption, Fig. 6 shows the excitation voltage required to reach 1 µm of vibration amplitude PtP.

Fig. 6: Experimental comparison between voltage of the controlled vibration system, and the response in the case of the double controlled system against temperature variation and vibration amplitude damping.

In this figure, we can verify that the voltage supply is higher when the system is not operated at the resonance. The difference between the excitation voltages is about 30% which indicates the importance of the two controllers when implemented in the tactile stimulator.

Friction control results

The previous control loops ensure the accuracy of the tactile stimulator itself on the point of view of physical variables such as vibration amplitude or energy consumption. However, the user has to be taken into account also on a tribological aspect. In this section, the interest is putted in the validation of the friction control strategy and its capability to impose a desired friction coefficient (COF) under the finger. The experiment consists on giving two references to the controller depending on the finger position. The difference between the two references in COF is imposed to be 0.1 in order to validate the controller robustness. Four subjects have been participated to this experiment, they were chosen according to their initial friction coefficient, which is different from a person to another, in order to test different operating points of the controller. The initial friction coefficient for the four subjects (ie without any vibration) is respectively between: 0.3 to 0.4, 0.4 to 0.5, 0.5 to 0.6 and 0.6 to 0.7. The COF reference given to the controller depending on each subject is receptively: 0.2 to 0.3, 0.3 to 0.4, 0.4 to 0.5 and 0.5 to 0.6. Fig. 4 illustrates the controlled friction coefficient between a finger and the plate for each subject. The desired friction coefficient is between the two values as explained in the foregoing depending on the finger position. Fig. 7 illustrates the experimental results for all the subjects, the x label represents the time in ms when the finger is altering between the two parts of the plate.

Fig. 7: Closed loop control of the friction coefficient as a function of time for the four participants. The dotted curves illustrate at least four try of touching the surface and the black continuous curve represents the average of all measurements for each subject.

When comparing the mean curve of the controlled friction for the 4 subjects, we find that the control strategy works well in static regime when imposing two references to the controller. In dynamic regime, there is a little difference when comparing the controlled friction for subjects; Tab. 1 shows a comparison in the dynamic regime of the rise time and the overshoot for each subjects. We find that the rise time is equal to 3.84 1.27 and the overshoot of 37.5 10.2 %.

Conclusion

In this paper, we proposed a new high fidelity smart tactile stimulator called "SMARTTAC". An approach to control the wave amplitude of the vibration was presented. The objective is to make the effect of the tactile feedback phenomenon insensitive to external disturbances. For that aim, a new modelling has been developed analytically due to an accurate choice of the d-q frame. Identification has been performed and then, a simultaneous double closed loop control has been implemented. The tracking of the resonant shows a high robustness against the temperature variation. The capability of this tactile stimulator to control in real time the friction coefficient was also exposed by proposing a direct control of the friction thanks to the sensors implementation in the stimulator. Tribological validation of the COF control was proved by repeated experiment with four subjects showing a stable and desired level of friction with a low response time of about 4 ms.

  Fig 4 illustrates the operated closed loop control strategy. High fidelity closed loop controlled friction in SMARTTAC tactile stimulator BEN MESSAOUD Wael EPE'15 ECCE Europe ISBN: 9789075815238 and CFP15850-USB P.5

Fig. 4 :

 4 Fig. 4: The proposed control stra loop which is disturbed by the va normal force, its velocity and dis loop acts on the desired vibration

  Control of vibration and frequency High fidelity closed loop controlled friction in SMARTTAC tactile stimulator BEN MESSAOUD Wael EPE'15 ECCE Europe ISBN: 9789075815238 and CFP15850-USB P.7

  High fidelity closed loop controlled friction in SMARTTAC tactile stimulator BEN MESSAOUD Wael EPE'15 ECCE Europe ISBN: 9789075815238 and CFP15850-USB P.8

Table I : The difference between the rise time and the overshoot between the four subjects.

 I

		Subject 1	Subject 2	Subject 3	Subject 4	Mean
	Rise time (ms)	3.37	3.41	2.87	5.72	3.84
	Overshoot (%)	27.6	44.4	48	30	37.5

© assigned jointly to the European Power Electronics and Drives Association & the Institute of Electrical and Electronics Engineers (IEEE)

Acknowledgements

The authors gratefully acknowledge financial support from the Nord Pas de Calais and Alsace regions. Special thanks to IRCICA and INRIA institutions for supporting this work.


Guillaume Burel

Gilles Dowek

Ying Jiang

A completion method to decide reachability in rewrite systems

The Knuth-Bendix method takes in argument a finite set of equations and rewrite rules and, when it succeeds, returns an algorithm to decide if a term is equivalent to another modulo these equations and rules. In this paper, we design a similar method that takes in argument a finite set of rewrite rules and, when it succeeds, returns an algorithm to decide not equivalence but reachability modulo these rules, that is if a term reduces to another. As an application, we give new proofs of the decidability of reachability in finite ground rewrite systems and in pushdown systems.

INTRODUCTION

The Knuth-Bendix method [START_REF] Knuth | Simple word problems in universal algebras[END_REF][START_REF] Huet | Confluent reductions: abstract properties and applications to term rewriting systems[END_REF] takes in argument a finite set of equations and rewrite rules and, when it succeeds, returns an algorithm to decide if a term is equivalent to another modulo these equations and rules. In this paper, we design a similar method that takes in argument a finite set of rewrite rules and, when it succeeds, returns an algorithm to decide not equivalence but reachability modulo these rules, that is if a term reduces to another.

As an application, we give new proofs of the decidability of reachability in finite ground rewrite systems [START_REF] Dauchet | Decidability of confluence for ground term rewriting systems[END_REF] and in pushdown systems [START_REF] Bouajjani | Reachability analysis of pushdown automata: Application to modelchecking[END_REF].

Like the Knuth-Bendix method, this method proceeds by completing a finite rewrite system into an equivalent one, by adding derivable rules. In the completed system, when a proposition t -→ * u has a proof, it also has a proof of the form t = t 0 -→ t 1 ... -→ t p = w = u q -→ ... u 1 -→ u 0 = u where t 0 t 1 ... t p and u q ≺ ... ≺ u 1 ≺ u 0 for some reduction order ≺, that is a proof formed with a decreasing sequence followed by an increasing one. We can write such a proof

t = t 0 u 0 = u t 1 - u 1
t p = w = u q -using the unusual convention to write the larger terms for the order ≺ on the top of the diagram and the smaller ones on the bottom, hence drawing an arrow oriented from the bottom to the top, when a smaller term reduces to a larger one.

In order to transform proofs into proofs of this form, we should consider critical steps of the form u t v with t ≺ u v and add a rule reducing directly t to v, avoiding the detour via u. If the reduction from t to u uses a rule l 1 -→ r 1 and that from u to v a rule l 2 -→ r 2 , the terms r 1 and l 2 would have to be compared, to determine if one term unifies with a subterm of the other. We would therefore need to design a forward completion method that compares the lefthand side of a rule with the right-hand side of another.

An alternative method is to reverse the rules whose left-hand side is smaller than the right-hand side, and keep track that such reversed rules must be used backwards. Thus, we distinguish two kinds of rules: negative rules that are as usual, and positive rules that must be used backwards: u 1 -→ + u 2 means that u 2 -→ u 1 in the original system and u 2 ≺ u 1 . This way, in the completed system, when a proposition t -→ * u has a proof, it also has a proof of the form

t = t 0 u 0 = u t 1 - - u 1 + t p = w = u q + --
The Computer Journal, Vol. ??, No. ??, ???? G. Burel, G. Dowek, and Y. Jiang and the critical steps have the form

u t + v - -
so that only left-hand side of rules need to be compared. In the completed system, the proposition t -→ * u has a proof if and only if there exists a term w such that t -→ * -w and u -→ * + w. Thus reachability boils down to the existence of a common reduct. In a terminating system, reachability is obviously decidable and easy to check because reduction trees are always finite. In the same way, in a terminating system, the existence of a common reduct of two terms is decidable and easy to check because reduction trees are finite.

The reader familiar with polarized sequent calculus modulo theory [START_REF] Dowek | What is a theory?[END_REF][START_REF] Dowek | Polarized resolution modulo[END_REF] will remark that many ideas in this paper, in particular the idea to distinguish two kinds of rules, come from this calculus. But the paper is presented independently of polarized sequent calculus modulo theory.

POLARIZED REWRITE SYSTEMS

We consider a finite set constants a, b, ... a finite set of function symbols f , g, ... and an infinite set of variables, x, y, X, Y , ... Terms, substitutions, rewrite rules, and rewrite systems are defined as usual.

A rewrite rule l -→ r is said to be left-linear if the term l is linear in each of its variables, that is if each variable of l occurs exactly once in l.

A context C[X 1 , ..., X n ] is an ordered pair formed with a term C, and a sequence of variables X 1 , ..., X n each occurring exactly once in C.

The term

(t 1 /X 1 , ..., t n /X n )C is written C[t 1 , ..., t n ].
Definition 2.1 (Polarized rewrite system). A polarized rewrite system P is a pair P -, P + of rewrite systems. The rules of P -are called negative and are written l -→ -r, the rules of P + are called positive and are written l -→ + r. The one step reduction relations -→ -and -→ + are defined as usual: t -→ -u (resp. t -→ + u) if there exists a negative rule l -→ -r (resp. a positive rule l -→ + r), a context C[X] and a substitution σ,

such that t = C[σl] and u = C[σr]. Definition 2.2 (The relation -→). The relation -→ is -→ -∪ + ←-, that is t -→ u if t -→ -u or u -→ + t.
Definition 2.3 (Proof). Let P be a polarized rewrite system. A proof (or a reduction sequence) in P is a sequence of terms t 0 , t 1 , ..., t n , such that for all i, t i -→ t i+1 , that is t i -→ -t i+1 or t i+1 -→ + t i .

A proof is a proof of the proposition t -→ * u if t = t 0 and u = t n . Definition 2.4 (Polarization). The polarized rewrite system P is said to be a polarization of a non-polarized rewrite system R if

• for each rule l -→ r of R, the system P contains either the rule l -→ -r or the rule r -→ + l, • for each rule l -→ -r of P, the system R contains the rule l -→ r, • and, for each rule l -→ + r of P, the system R contains the rule r -→ l.

Proposition 2.1. Let R be a rewrite system, and P be a polarization of R, then a proposition t -→ * u has a proof in R if and only if it has a proof in P.

Example 1. Let R be the rewrite system

f (g(a)) -→ b f (f (a)) -→ c b -→ g(g(f (f (a)))) c -→ h(a)
The polarized rewrite system P

f (g(a)) -→ -b f (f (a)) -→ -c g(g(f (f (a)))) -→ + b h(a) -→ + c is a polarization of R. The proposition f (g(a)) -→ * g(g(h(a))) has the following proof in R f (g(a)) g(g(f (f (a)))) g(g(h(a))) b - - g(g(c))
-and the following proof in P

f (g(a)) g(g(f (f (a)))) g(g(h(a))) b + - - g(g(c)) + -- Definition 2.5 (Termination). A polarized rewrite system is terminating if the relation -→ -∪ -→ + is well-founded.

Note that this does not imply that the relation

-→, that is -→ -∪ + ←-, is well-founded.
Definition 2.6 (Reduction order). A reduction order ≺ is an order such that

• if t ≺ u, then for all function symbols f and terms t 1 , ..., t i-1 , t i+1 , ..., t n f (t 1 , ..., t i-1 , t, t i+1 , ..., t n ) ≺ f (t 1 , ..., t i-1 , u, t i+1 , ..., t n ) • if t ≺ u, then for all substitutions σ σt ≺ σu • ≺ is well-founded.
Proposition 2.2. Let ≺ be a reduction order and P be a polarized rewrite system such that l r for each rule l -→ -r or l -→ + r of P. Then, the system P terminates.

CUT-ELIMINATION

Definition 3.1 (Cut). A cut (or a peak) in a proof π = t 0 , t 1 , ..., t n is a sub-sequence t i-1 , t i , t i+1 such that t i t i-1 + t i+1 - -
A proof is cut-free (or a valley proof) if it contains no cuts, that is if it is formed with a sequence of negative steps followed by a sequence of positive steps. A polarized rewrite system has the cut-elimination property (or is confluent) if every proposition t -→ * u that has a proof has a cut-free proof.

Example 2. In the polarized rewrite system P of Example 1, the proposition f (g(a)) -→ * g(g(h(a))) has the proof

f (g(a)) g(g(f (f (a)))) g(g(h(a))) b + - - g(g(c)) + -- but no cut-free proof.
The reader familiar with polarized sequent calculus modulo theory will remark that the proposition t -→ * u has a proof (resp. a cut-free proof) in P if and only if the sequent P (t) P (u), where P is a predicate symbol, has a proof (resp. a cut-free proof) in polarized sequent calculus modulo P. Proposition 3.1. Let P be a terminating finite polarized rewrite system. Then, the existence of a cutfree proof in P of a proposition t -→ * u is decidable.

Proof. The proposition t -→ * u has a cut-free proof if and only if the reducts of t in P -and those of u in P + have a term in common. As P terminates, both reduction trees are finite.

Definition 3.2 (Proof reduction).

A proof π reduces to a proof π , if π is obtained by replacing a cut in π by a cut-free proof, that is if

π = t 0 , ..., t i-1 , t i , t i+1 , ..., t n with t i t i-1 + t i+1 - - and π = t 0 , ..., t i-1 = u 0 , u 1 , ..., u p = w = v q , ..., v 1 , v 0 = t i+1 , ..., t n with t i-1 = u 0 t i+1 = v 0 u 1 -- v 1 + u p = w = v q + -- Definition 3.3 (Local confluence). A polarized rewrite system is locally confluent if every cut is reducible, that is if for each proof t u + v - -
there exists a proof

u = u 0 v = v 0 u 1 - - v 1 + u p = w = v q + --
Newman's lemma can be seen as a termination lemma for proof-reduction [START_REF] Dowek | Confluence as a cut elimination property[END_REF], following an idea that is already in [START_REF] Jouannaud | Completion of a set of rules modulo a set of equations[END_REF]. Proposition 3.2 (Newman). If P is terminating and locally confluent, then it has the cut-elimination property.

Proof. As P is terminating, the transitive closure of the relation -→ - ∪ -→ + is a well-founded order. Thus, its multiset extension < [START_REF] Dershowitz | Proving termination with multiset orderings[END_REF] is also well-founded. A proof-reduction step replaces the multiset {t 1 , ..., t i-1 , t i , t i+1 , ..., t n } with the multiset {t 1 , ..., t i-1 , u 1 , ..., u p-1 , w, v q-1 , ..., v 1 , t i+1 , ..., t n } and {t 1 , ..., t i-1 , u 1 , ..., u p-1 , w, v q-1 , ..., v 1 , t i+1 , ..., t n } < {t 1 , ..., t i-1 , t i , t i+1 , ..., t n } because each term u 1 , ..., u p-1 , w, v q-1 , ..., v 1 is smaller than t i . Thus, proof-reduction terminates.

Finally, as P is locally confluent, an irreducible proof contains no cuts. Definition 3.4 (Critical pair). A critical pair is a pair of terms of the form

• σr 1 , (σC)[σr 2 ] , where l 1 -→ -r 1 is a negative rule, l 2 -→ + r 2 is a positive rule, C[X] is a context, l 1 is
a term that is not a variable, and σ is a substitution, such that X does not occur in σ, l 1 = C[l 1 ], and σ is the most general unifier of l 1 and l 2 ,

• or (σC)[σr 1 ], σr 2 , where l 1 -→ -r 1 is a negative rule, l 2 -→ + r 2 is a positive rule, C[X] is a context, l 2 is a term that is not a variable, and

The Computer Journal, Vol. ??, No. ??, ???? G. Burel, G. Dowek, and Y. Jiang σ is a substitution, such that X does not occur in σ, l 2 = C[l 2 ], and σ is the most general unifier of l 1 and l 2 .

A critical pair u, v is joinable if there exists a term w, such that the propositions u -→ * + w and v -→ * -w are provable.

We now would like to prove that if all the critical pairs of a polarized rewrite system P are joinable, then P is locally confluent. Unfortunately, this property does not hold in general, as shown by the following counterexample

f (x, x) -→ -g(x) a -→ + b
that contains no critical pairs, but that is not locally confluent, as the cut

f (a, a) f (a, b) + g(a)
-cannot be reduced: the term g(a) reduces positively to g(b) only, and the term f (a, b) cannot be reduced negatively.

It indeed reduces to f (b, b), but not negatively.

N. Hirokawa [START_REF] Hirokawa | Commutation and signature extension[END_REF] has found a similar counter-example independently, in a different context. Fortunately, this property holds for left-linear rewrite systems.

We start by recalling two well-known classification lemmas [START_REF] Knuth | Simple word problems in universal algebras[END_REF][START_REF] Huet | Confluent reductions: abstract properties and applications to term rewriting systems[END_REF].

Proposition 3.3. Let C 1 [X] and C 2 [
Y ] be contexts, and u 1 and u 2 be terms such that

C 1 [u 1 ] = C 2 [u 2 ] then • either the occurrences of X and Y are disjoint, that is there exists a context D[X, Y ] such that C 1 [X] = D[X, u 2 ] and C 2 [Y ] = D[u 1 , Y ] C 1 C 2 D
• or the occurrence of X is higher than that of Y , that is there exists

a context D[Y ] such that C 2 [Y ] = C 1 [D[Y ]] C 1 C 2 D C 1
• or the occurrence of Y is higher than that of X, that is there exists a

context D[X] such that C 1 [X] = C 2 [D[X]] D C C C 1 2 2
Proposition 3.4. Let t and u be terms, τ be a substitution and D[Y ] be a context such that the variable Y does not occur in τ and τ t = D[u]. Then

• either the occurrence of Y in D is not an occurrence of t, that is there exist a variable x and contexts

E 1 [X] and E 2 [Y ], such that t = E 1 [x], τ x = E 2 [u], and D[Y ] = (τ E 1 )[E 2 [Y ]] τ 1 2 t D E E
• or the occurrence of Y in D is an occurrence of t, that is there exist a context E and a term t such that t = E[t ] and

D[Y ] = (τ E)[Y ] τ t t' D E
Proposition 3.5. If all the critical pairs of a leftlinear polarized rewrite system P are joinable, then P is locally confluent.

Proof. Consider three terms t, u, and v such that

t u + v - -
where t reduces to u by a rule l 1 -→ + r 1 , and to v by a rule l 2 -→ -r 2 . As variables in rules may be renamed we can assume, without loss of generality, that l 1 and l 2 have no variables in common.

There exist two contexts C 1 [X] and C 2 [Y ] and a substitution τ , such that X and Y do not occur in

τ , t = C 1 [τ l 1 ] = C 2 [τ l 2 ], u = C 1 [τ r 1 ], and v = C 2 [τ r 2 ]
. Thus, by Proposition 3.3, either there exists

a context D[X, Y ] such that C 1 [X] = D[X, τ l 2 ] and C 2 [Y ] = D[τ l 1 , Y ], or there exists a context D[Y ] such that C 2 [Y ] = C 1 [D[Y ]], or there exists a context D[X] such that C 1 [X] = C 2 [D[X]
]. We consider these three cases.

• If there exists a context D[X, Y ] such that C 1 [X] = D[X, τ l 2 ] and C 2 [Y ] = D[τ l 1 , Y ] C 1 C 2 D we have u = D[τ r 1 , τ l 2 ], and v = D[τ l 1 , τ r 2 ], let w = D[τ r 1 , τ r 2 ]. We have u v w * + * - - • If there exists a context D[Y ] such that C 2 [Y ] = C 1 [D[Y ]] C 1 C 2 D C 1 we have t = C 1 [τ l 1 ] = C 1 [D[τ l 2 ]], u = C 1 [τ r 1 ], and v = C 1 [D[τ r 2 ]]. As C 1 [τ l 1 ] = C 1 [D[τ l 2 ]], we have τ l 1 = D[τ l 2 ]
. Therefore, by Proposition 3.4, either there exist a variable x and contexts

E 1 [X] and E 2 [Y ], such that l 1 = E 1 [x], τ x = E 2 [τ l 2 ],
and

D[Y ] = (τ E 1 )[E 2 [Y ]],
or there exist a context E[Y ] and a term l 1 such that l 1 = E[l 1 ] and

D[Y ] = (τ E)[Y ].
We consider these two cases.

-

If l 1 = E 1 [x], τ x = E 2 [τ l 2 ],
and

D[Y ] = (τ E 1 )[E 2 [Y ]
], then we let τ = τ |V\{x} , and we have τ = (τ , E 2 [τ l 2 ]/x). The term l 1 is linear and

x does not occur in E 1 , thus τ E 1 = τ E 1 . Let w = C 1 [(τ , E 2 [τ r 2 ]/x)r 1 ]. We have u = C 1 [τ r 1 ] = C 1 [(τ , E 2 [τ l 2 ]/x)r 1 ] -→ * -C 1 [(τ , E 2 [τ r 2 ]/x)r 1 ] = w and v = C 1 [D[τ r 2 ]] = C 1 [(τ E 1 )[E 2 [τ r 2 ]]] = C 1 [(τ E 1 )[E 2 [τ r 2 ]]] = C 1 [(τ , E 2 [τ r 2 ]/x)(E 1 [x])] = C 1 [(τ , E 2 [τ r 2 ]/x)l 1 ] -→ * + C 1 [(τ , E 2 [τ r 2 ]/x)r 1 ] = w u v w * + * - - - If l 1 = E[l 1 ] and D[Y ] = (τ E)[Y ], then τ l 1 = (τ E)[τ l 1 ] = D[τ l 1 ]. As we have τ l 1 = D[τ l 2 ], we get D[τ l 1 ] = D[τ l 2 ], thus τ l 1 = τ l 2 .
Let σ be the most general unifier of l 1 and l 2 and η such that τ = η • σ.

We have

u = C 1 [τ r 1 ] = C 1 [ησr 1 ] and v = C 1 [(ησE)[ησr 2 ]] = C 1 [η((σE)[σr 2 ])].
We know that the critical pair σr 1 , (σE)[σr 2 ] closes on a term, say

w 0 . Let w = C 1 [ηw 0 ]. We have u v w * + * - -
• The third case is similar to the second.

Definition 3.5 (Polarized Knuth-Bendix method). Let P be a left-linear finite polarized rewrite system and ≺ a reduction order, such that l r for each rule l -→ -r or l -→ + r of P.

While there is a non-joinable critical pair

t u + v - -
if u v and u -→ -v is a left-linear rewrite rule, add this rule to close the critical pair, if v u and v -→ + u is a left-linear rewrite rule, add this rule to close the critical pair, otherwise fail.

Proposition 3.6. Let P be a left-linear finite polarized rewrite system and ≺ a reduction order, such that l r for each rule l -→ -r or l -→ + r of P. If the polarized Knuth-Bendix method applied to P succeeds, then reachability in P is decidable.

Proof. Let P be the left-linear finite polarized rewrite system built by the polarized Knuth-Bendix method. The rules of P are all derivable in P, thus a proposition t -→ * u has a proof in P if and only if it has a proof in P .

As all the critical pairs of P are joinable, by Proposition 3.5, P is locally confluent. By construction, l r for each rule l -→ -r or l -→ + r of P . Thus, by Proposition 2.2, P terminates. By Proposition 3.2, as P is locally confluent and terminating, it has the cut-elimination property.

Thus, a proposition t -→ * u has a proof in P if and only if it has a proof in P if and only if it has a cutfree proof in P . And, by Proposition 3.1, the existence of a cut-free proof for a proposition t -→ * u in P is decidable.

Example 3. Let P be the system defined in Example 1 and ≺ be the Knuth-Bendix order [START_REF] Knuth | Simple word problems in universal algebras[END_REF] with an equal weight 1 for all symbols and any precedence. For all rules l -→ -r or l -→ + r of P, we have l r. The only non-joinable critical pair is Let R be a (non-polarized) finite rewrite system, P be a polarization of R and ≺ be a reduction order. If P is a left-linear polarized rewrite system, for all rules l -→ -r or l -→ + r of P, l r, and the polarized Knuth-Bendix method applied to P succeeds, then reachability in R is decidable.

Proof. From Propositions 2.1 and 3.6.

The Computer Journal, Vol. ??, No. ??, ???? G. Burel, G. Dowek, and Y. Jiang

GROUND FINITE REWRITE SYSTEMS

A ground rewrite system is a rewrite system such that for all rules l -→ r, both terms l and r are ground. D. Lankford [START_REF] Lankford | Canonical inference[END_REF][START_REF] Dershowitz | Completion and its applications[END_REF] has observed that if ≺ is the Knuth-Bendix order with an equal weight 1 for all symbols and any precedence, R is a finite ground rewrite system, R is the equivalent system obtained by removing the rules of the form l -→ l and reversing the rules l -→ r such that l ≺ r into r -→ l, then the Knuth-Bendix method always succeeds on R , and therefore equivalence in R is decidable.

We now want to prove that, in a similar way, reachability in a finite ground rewrite system is decidable [START_REF] Dauchet | Decidability of confluence for ground term rewriting systems[END_REF].

Theorem 4.1 (Dauchet-Tison). Let R be a finite ground rewrite system. Then, the existence of a proof in R of a proposition t -→ * u is decidable.

Proof. Let ≺ be the Knuth-Bendix order with an equal weight 1 for all symbols and any precedence. This order is a reduction order and it is total on ground terms.

Without loss of generality, we can assume that R does not contain trivial rules of the form l -→ l. Let P be the polarization of R obtained by transforming each rule l -→ r of R such that l r, into a negative rule l -→ -r, and reversing each rule l -→ r such that l ≺ r into a positive rule r -→ + l. By construction, l r for each rule l -→ -r or l -→ + r of P.

Let T be the finite set containing the left-hand sides of the rules of P and T be the set of ground terms t such that there exists a term u in T such that t ≺ u or = u. As, for the Knuth-Bendix order, if u is a ground term, the set of ground terms t such that t ≺ u is always finite, the set T is finite.

Then, the polarized Knuth-Bendix method applied to P generates rules whose left-hand sides and right-hand sides are in T . As there is only a finite number of such rules, the polarized Knuth-Bendix method applied to P terminates successfully.

Note that the original proof based on the construction of automata recognizing left-hand sides and right-hand sides of rules [START_REF] Dauchet | Decidability of confluence for ground term rewriting systems[END_REF]

f (s 1 ) -→ s s -→ s s 1 -→ b s -→ g(s 1 )
where

f (a) -→ f (s 1 ) -→ s -→ s -→ g(s 1 ) -→ g(b)
and then reversing the two rules

s 1 -→ b s -→ g(s 1 )
The first step is in fact not needed.

PUSHDOWN SYSTEMS

As another corollary of our result, we also get the decidability of reachability for pushdown systems [START_REF] Bouajjani | Reachability analysis of pushdown automata: Application to modelchecking[END_REF].

Definition 5.1 (Pushdown system). Consider a language containing a set S of unary function symbols called stack symbols, a set Q of unary function symbols called states and a constant ε. A pushdown system is a finite rewrite system with rules of the form: pop rules p(γ(x)) -→ q(x) where γ is a stack symbol and p and q are states, push rules p(x) -→ q(γ(x))

where γ is a stack symbol and p and q are states, and neutral rules p(x) -→ q(x)

where p and q are states.

Theorem 5.1 (Bouajjani-Esparza-Maler). Let R be a pushdown system. Then, the existence of a proof in R of a proposition t -→ * u is decidable.

Proof. Consider a total precedence on function symbols such that stack symbols are larger than states and let ≺ be the lexicographic path order [START_REF] Kamin | Attempts for generalizing the recursive path ordering[END_REF] relative to this precedence. The order ≺ is a reduction order and p(γ(x)) q(x) for all γ, p and q.

Pop rules are polarized as p(γ(x)) -→ -q(x) push rules as q(γ(x)) -→ + p(x)

and neutral rules according to the precedence. Critical pairs have the form p(γ(x)) q(x) + r(x)

--p(γ(x)) q(γ(x)) + r(x)

-p(γ(x)) q(x) + r(γ(x))

- - p(x) q(x) + r(x)
-all these critical pairs are closed by adding a pop rule, a push rule, or a neutral rule. As there are a finite number of such rules, the polarized Knuth-Bendix method always terminates successfully. Thus reachability in pushdown systems is decidable.

Note that we get, in this way, the decidability of reachability in pushdown systems, but not of alternating pushdown systems [START_REF] Bouajjani | Reachability analysis of pushdown automata: Application to modelchecking[END_REF], that requires moving from polarized rewrite systems to polarized sequent calculus modulo [START_REF] Burel | Automata, resolution and cut-elimination[END_REF].

  with the rule g(g(c)) -→ + bLet P be the system obtained by adding this rule to P. The proposition f (g(a)) -→ * g(g(h(a))) has the proof in P

  also uses implicitly the idea of reversing rewrite rules. For instance, with the rule f (a) -→ g(b), it builds an automaton recognizing f (a) in s a -→ s 1 f (s 1 ) -→ s another recognizing g(b) in s b -→ s 1 g(s 1 ) -→ s and takes the rewrite rule s -→ s This construction can be decomposed in two steps, one transforming the rule f (a) -→ g(b) into the rewrite system a -→ s 1

ACKNOWLEDGEMENTS

The authors want to thank Nao Hirokawa, Gérard Huet, Jean-Pierre Jouannaud, Claude Kirchner, and Vincent van Oostrom, for useful discussions on this paper and bibliographical indications. This work is supported by the ANR-NSFC project LOCALI (NSFC 61161130530 and ANR 11 IS02 002 01) and the Chinese National Basic Research Program (973) Grant No. 2014CB340302.