L Krapivsky

Kirone Mallick

Chikashi Arita

P L Krapivsky

Variational calculation of transport coefficients in diffusive lattice gases

published or not. The documents may come

Variational calculation of transport coefficients in diffusive lattice gases

Chikashi Arita, P.

I. INTRODUCTION

In systems composed of a huge number of interacting particles, a continuous macroscopic behavior emerges after space and time variables are suitably rescaled. In this hydrodynamic limit, the intrinsic granularity of the basic constituents is lost and the system is described by continuous fields (such as matter or charge densities, currents, magnetization, etc.), coupled by partial differential equations. The program of developing a mathematical theory of hydrodynamic limits for general systems was posed by Hilbert [START_REF]Boltzmann's work on the principles of mechanics suggests the problem of developing mathematically the limiting processes, there merely indicated, which lead from the atomistic view to[END_REF] in his sixth problem, and is still far from completion.

Over the past 30 years, significant progress in deriving hydrodynamic limits has been achieved in the realm of lattice gases with stochastic microscopic dynamics [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF][START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF][START_REF] Liggett | Stochastic Interacting Systems: Contact, Voter, and Exclusion Processes[END_REF][START_REF] Katz | [END_REF][6][7][START_REF] Schmittmann | Statistical Mechanics of Driven Diffusive Systems[END_REF][START_REF] Schütz | Exactly Solvable Models for Many-Body Systems Far From Equilibrium[END_REF][START_REF] Blythe | [END_REF][START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF][START_REF] Derrida | [END_REF][13][14][15]. In particular, for stochastic lattice gases with simple conservative interactions such as exclusion processes it was shown [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF][START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF] that a coarse-grained density ρ(x, t) satisfies the macroscopic conservation equation, ∂ t ρ = -∂ x J, with local current J(x, t) given by Fick's law, J = -D(ρ)∂ x ρ. The diffusivity D(ρ) depends on the microscopic dynamical rules and calculating, it is a very challenging problem.

The hydrodynamic limit describes the deterministic evolution of the density field ρ(x, t) defined as a local, coarse-grained empirical average over microscopic configurations. This is analogous to the law of large numbers [6]. The next step is to investigate fluctuations around the average, i.e., to find a property analogous to the central limit theorem for interacting lattice gases. A non-rigorous but physically well-motivated approach is to consider the density and current as stochastic fields coupled by mass conservation, ∂ t ρ = -∂ x J, and to add a random contribution to the constitutive equation that relates current to density: J = -D(ρ)∂ x ρ + σ(ρ)ξ(x, t) .

(1)

Here ξ(x, t) is a Gaussian white noise. The amplitude of the noise depends on the second transport coefficient σ(ρ) known as conductivity (or mobility). Similarly to the diffusivity, the conductivity is a function of the local density. The conductivity σ(ρ) depends on the microscopic rules. The two transport coefficients D(ρ) and σ(ρ) are difficult to calculate. If one of them is known, however, it is usually simple to determine the other due to the Einstein relation. For lattice gases close to equilibrium, the Einstein relation acquires a simple form

2D(ρ) σ(ρ) = d 2 F(ρ) dρ 2 , (2)
where F(ρ) is the equilibrium free energy density per site. Although a rigorous derivation of the stochastic partial differential equation obtained from [START_REF]Boltzmann's work on the principles of mechanics suggests the problem of developing mathematically the limiting processes, there merely indicated, which lead from the atomistic view to[END_REF] is lacking, the predictions agree with available exact results established for special lattice gases, see e.g. [START_REF] Derrida | [END_REF]13]. There are also independent mathematical arguments supporting the large deviation principle implied by the stochastic Langevin equation based on [START_REF]Boltzmann's work on the principles of mechanics suggests the problem of developing mathematically the limiting processes, there merely indicated, which lead from the atomistic view to[END_REF]. This is the starting point of the macroscopic fluctuation theory that describes diffusive interacting particle models in infinite domains and in finite systems connected to reservoirs at different temperatures (or chemical potentials) [14,15]. As long as the local equilibrium is satisfied, the macroscopic fluctuation theory is applicable to far-from-equilibrium regimes.

Thus, an understanding of macroscopic behaviors, both the deterministic (hydrodynamic) part and fluctuations around it, requires knowledge of two transport coefficients: D(ρ) and σ(ρ). The Einstein relation [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF] implies that it suffices to determine one coefficient. We focus on D(ρ) which is especially important since it governs the hydrodynamic behavior.

In recent years, a significant effort has been devoted to the calculation of transport coefficients for various lattice gases. When a stochastic lattice gas satisfies a special property known as the gradient condition [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF][START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF], the computations become feasible. The gradient property states that the microscopic current is the gradient of a local function (i.e., loosely speaking, Fick's law is already valid at the microscopic level). The simplest lattice gas obeying the gradient property is a collection of non-interacting random walkers-in this case, D = 1 and σ = 2ρ. The simplest interacting gradient lattice gas is the symmetric simple exclusion process [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF]16] for which D = 1 and σ = 2ρ(1 -ρ). Other gradient lattice gases for which the diffusivity has been computed include the Katz-Lebowitz-Spohn model with symmetric hopping [START_REF] Katz | [END_REF]17,18], repulsion processes [19], a lattice gas of leap-frogging particles [20,21], and an exclusion process with avalanches [22]. In these models, an exact expression for the diffusivity can also be derived by a "perturbation approach": one writes the current at the discrete lattice level and performs a continuous limit assuming the density field to be slowly varying.

Generic interacting particle processes do not satisfy the gradient condition; the computation of the diffusivity in such gases appears intractable. Nevertheless, there exists an exact variational formula for the diffusivity D(ρ) derived by Varadhan and Spohn [16] (see also [START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF]23,24]). This rather abstract formula is valid for general lattice gases regardless of the gradient condition. It expresses D(ρ) as a minimum of a functional over certain classes of functions. In simple cases, the functional is quadratic and the minimization gives a set of linear equations; generally it is unclear how to determine the minimum because the function space is infinite-dimensional.

In this work we demonstrate that the abstract variational formula for the diffusivity can be used as a tool to derive explicit (albeit approximate) formulas. We implement a systematic approximation procedure for D(ρ) by minimizing over finite-dimensional subspaces of the infinite-dimensional function space. The simplest approximation gives exact results for gradient lattice gases. For general lattice gas, this iterative scheme can be carried out analytically as far as one wishes, although the complexity of calculations increases rapidly with the dimensionality of the subspace. The precision greatly improves after each step.

The general idea of approximately solving a variational problem is widely known in science and engineering. For instance, it underlies the Ritz method. In quantum mechanics, it is known as the variational method used, e.g., in finding approximations to the ground state, and also to excited states. This venerable idea has not yet been applied to the Varadhan-Spohn formula for the diffusivity, mostly because it is little known, and it has the reputation of being a very abstract object that makes concrete calculations difficult. Furthermore, an approximate variational procedure based on the Varadhan-Spohn formula requires very long calculations even for the simplest nongradient lattice gases. Therefore one would like to choose a nongradient lattice gas which is natural and sufficiently simple to be amenable to analysis. Exclusion processes (lattice gases with at most one particle per site) appear to be good candidates. They are widely used in conceptual developments such as testing far-from-equilibrium behaviors (see [START_REF] Derrida | [END_REF][13][14][15][16]), and in various applications (see [25,26] FIG. 1. Illustration of the GEP in one dimension, where each site is occupied by at most two particles. The arrows indicate possible transitions with rates (3). Two particles occupying the same site jump independently, each with a rate 1 2 p2s.

and references therein). The basic example, namely the symmetric simple exclusion process, is gradient. In nongradient exclusion processes, the range of hopping is increased or there are interactions between particles occupying neighboring sites. In such situations the dimensionality of the subspace also increases-this makes computations unwieldy, and even more so if equilibrium is not given by a product measure [START_REF]product measure" we mean that in the equilibrium correlation functions factorize[END_REF].

A rather simple nongradient lattice gas with product measure is the 2-GEP, a generalized exclusion process (GEP) in which each site can host at most two particles. The 2-GEP is additionally parametrized by hopping rates depending on the occupancy levels. The 2-GEP is the first member in the family of k-GEP, where k ≥ 2 is the maximal occupancy. The k-GEPs are nongradient for all k ≥ 2 and for generic hopping rates; the notable exception is the misanthrope process [START_REF] Cocozza-Thivent | [END_REF]29].

The GEPs have been investigated in a number of studies [START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF]23,24,[30][31][32]. We overview their basic properties in Sec. II. In Sec. III we present the Varadhan-Spohn formula for the k-GEPs and we develop an iterative procedure allowing us to find increasingly better upper bounds for diffusivity. For the 2-GEP, each iteration of the variational procedure improves the precision by an order of magnitude. In Sec. IV, we investigate the 2-GEP on a finite interval with open boundaries connected to reservoirs. We conclude in Sec. V.

II. GENERALIZED EXCLUSION PROCESSES

We now define generalized exclusion processes (GEPs) and recall some of their basic properties at equilibrium. We start with the general k-GEP and then discuss the 2-GEP which is the main focus of our study. In the following we consider lattice gases in one dimension and always assume that particles undergo nearest-neighbor symmetric hopping.

A. Definition of k-GEPs

For the k-GEP, each site can accommodate at most k particles (see Fig. 1), i.e. the local occupation variables satisfy τ i ∈ {0, 1, . . . , k}. Jumps between adjacent sites are allowed only to sites with fewer than k particles and the rates can depend on the occupancies of the sites: (τ i , τ i+1 ) = (r, s) → (r -1, s + 1) (rate p rs ), (τ i , τ i+1 ) = (r, s) → (r + 1, s -1) (rate p sr ).

(

) 3
The hopping is possible only from occupied sites and impossible into maximally occupied sites, so

p 0s = p rk = 0. ( 4
)
The one-site probabilities in the equilibrium state are conveniently expressed through the fugacity λ,

P[τ i = r] = W r = a r λ r Z , (5)
where Z is a normalization factor

Z = 0≤r≤k a r λ r . ( 6
)
Without loss of generality we can set a 0 = a 1 = 1, so that the fugacity is defined by the ratio

W 1 /W 0 = λ.
The coefficients a r depend on hopping rates (3). Note that the fugacity is implicitly determined by the density

ρ = τ i = λ d dλ ln Z . (7)
In many GEPs, equilibrium is characterized by a product measure. For such lattice gases, the compressibility is defined by χ = τ 2 i -ρ 2 , and it can be alternatively written as

χ = τ 2 i -ρ 2 = λ dρ dλ . (8)
For lattice gases with equilibrium having the product structure, the equilibrium free energy density per site also admits a simple general form (see e.g. [START_REF] Godrèche | Ageing and the Glass Transition[END_REF])

F = ρ ln λ -ln Z . (9)
Differentiating F with respect to ρ and using (7) we obtain dF dρ = ln λ. Differentiating again gives d 2 F dρ 2 = 1 χ , so that the Einstein relation (2) can be rewritten as

σ = 2χD. ( 10
) B. 2-GEP
For the 2-GEP, there are four independent generally non-vanishing hopping rates: p 10 , p 11 , p 20 , p 21 . The remaining five vanish: p 00 = p 01 = p 02 = p 12 = p 22 = 0.

Some special sets of rates have been considered in the literature, particularly

• Particle-uniform rates: Each particle hops with the same (unit) rate as long as the maximal occupancy constraint is obeyed (see e.g. [34]). Therefore the non-vanishing rates are

p 20 = p 21 = 2, p 10 = p 11 = 1. (11)
• Site-uniform rates: The occupancy of each site is updated with rate 1 as as long as the maximal occupancy constraint is obeyed (see e.g. [START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF]30,32])

p 20 = p 21 = p 10 = p 11 = 1. (12)
• The misanthrope process [START_REF] Cocozza-Thivent | [END_REF]:

p 20 = p 21 + p 10 . (13)
Whenever possible we shall consider the general 2-GEP with parameters p 10 , p 11 , p 20 , p 21 being arbitrary nonnegative numbers, but some results valid specifically for the rates ( 11), ( 12) or ( 13) will be emphasized. In particular, simulations have been performed for the 2-GEP with particle-uniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF].

We have three one-site probabilities in the equilibrium

W 0 = 1 Z , W 1 = λ Z , W 2 = a λ 2 Z ( 14
)
with a 2 = a for simplicity. The normalization factor (6) turns into a quadratic polynomial

Z = 1 + λ + aλ 2 (15)
for the 2-GEP. The product measure holds for the 2-GEP, viz. the probability of each configuration takes the product form [START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF] at equilibrium. In particular,

P[τ i = r, τ i+1 = s] = W r W s . (16)
The validity of the product measure can be seen by checking the detailed balance condition corresponding to the processes (r, s) ↔ (r -1, s + 1):

p rs W r W s = p r+1,s-1 W s-1 W r+1 . ( 17
)
For the 2-GEP, these equations are identities for most (r, s), the only exception is

p 20 W 2 W 0 = p 11 W 2 1 , (18)
which fixes the coefficient a in ( 14) to be

a = p 11 p 20 . (19)
The single-site probabilities W r , and generally the probabilities of arbitrary configurations which are the products of single-site probabilities, could have been functions of p 10 , p 11 , p 20 , p 21 . However, they actually depend only on the ratio (19).

For the 2-GEP, the density and the compressibility can be written as

ρ =W 1 + 2W 2 , (20)
χ =W 1 + 4W 2 -(W 1 + 2W 2 ) 2 . (21)
Combining ( 14) and (20) we determine an explicit formula for λ in terms of ρ λ = 2ρ

1 -ρ + 1 -(1 -4a)ρ(2 -ρ) (22)
and then we express the compressibility via density:

χ = ρ(2 -ρ) 1 -(1 -4a)ρ(2 -ρ) 1 + 1 -(1 -4a)ρ(2 -ρ) . (23)
We always assume that 0 < a < ∞. Peculiar behaviors may occur in the extreme cases of a = 0 and a = ∞. These subtleties are outlined in Appendix A. Possible qualitative differences between the 2-GEP and k-GEPs with k ≥ 3 are outlined in Appendix B.

III. BULK DIFFUSIVITY: A VARIATIONAL CALCULATION

A. Varadhan-Spohn formula for the diffusivity

We first review the Varadhan-Spohn variational formula giving the diffusion coefficient D(ρ) of a stochastic lattice gas. Generally, D(ρ) can be expressed, via the Green-Kubo formula, as an integral of a current-current correlation function [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF][START_REF] Katz | [END_REF]. The Green-Kubo expression can be rewritten as the solution of a variational problem (as shown by Spohn [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF] who attributes this variational approach to Varadhan). More precisely [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF]24]:

D = 1 2χ inf f Q(f ) . ( 24
)
Here χ is the compressibility and Q(f ) is a quadratic functional of the space of functions f which depend only on finite points of τ ("cylinder" functions). The expectation value • is taken with respect to the equilibrium measure on the configuration space, see Sec. II. The precise form of the functional Q(f ) depends on the microscopic dynamical rules of the process. To write down Q(f ) for the k-GEP, it is convenient to use some auxiliary notations. For any configuration τ = (. . . , τ -1 , τ 0 , τ 1 , τ 2 , . . .) we denote by τ 0→1 the configuration obtained from τ by making a single particle jump from site 0 to site 1 if such move move is permitted, otherwise τ 0→1 is identical to τ :

τ 0→1 = (. . . , τ 0 -1, τ 1 + 1, . . .) τ 0 ≥ 1, τ 1 < k τ otherwise . ( 25
)
When the particle hop is possible, two occupation numbers in τ 0→1 differ from the corresponding occupation numbers in τ , and only these numbers are explicitly shown in Eq. (25) . Similarly the configuration τ 0←1 is obtained from τ by making a particle hop from site 1 to site 0 if possible:

τ 0←1 = (. . . , τ 0 + 1, τ 1 -1, . . .), (26)
when τ 0 < k and τ 1 ≥ 1; otherwise τ 0←1 = τ . For any cylinder function f and j ∈ Z, we define f j by

f j (τ ) = f (T j τ ), (27)
where the translation operator T j shifts the configuration τ forward by j sites: (T j τ ) i = τ i-j . In other words, we shift the frame of f backward by j sites. For instance, if f depends only on the (τ

1 , τ 2 • • • , τ n ), that is, f (τ ) = f (τ 1 , τ 2 • • • , τ n ), (28)
the shifted function f j acts on the configuration space as

f j (τ ) = f (τ 1-j , τ 2-j • • • , τ n-j ). ( 29
)
Note that f 0 ≡ f . We can now define Q(f ) by its action on the space of cylinder functions f :

Q(f )(τ ) = p τ0τ1 1 - j∈Z f j (τ 0→1 ) -f j (τ ) 2 + p τ1τ0 -1 - j∈Z f j (τ 0←1 ) -f j (τ ) 2 . ( 30
)
The sum with respect to j in (30) contains a finite number of terms when the function f has a finite range. To appreciate this assertion consider a function of type (29). For j < 0 and j > n, the shifted function f j is not sensitive to the values of τ 0 and τ 1 , i.e., in Eq. ( 29) τ 0 and τ 1 do not appear. Therefore f j (τ 0→1 )-f j (τ ) = 0 for such j. Thanks to this property, one can replace j∈Z in (30) by 0≤j≤n . Only a finite number of variables, viz. (τ -n+1 , τ -n , • • • , τ n-1 , τ n ), appear in (30).

B. Iterative scheme

The Varadhan-Spohn formula ( 24) is a powerful theoretical tool, but it is unclear how to apply it even to simplest lattice gases such as the 2-GEP. One can try to obtain an upper bound for the diffusivity by restricting the space of functions. Let us consider the subspace F n of functions depending on the configuration of n adjacent sites (τ 1 , τ 2 , . . . , τ n ) and use a restricted version of the Varadhan-Spohn formula

D n = 1 2χ min f ∈Fn Q(f ) . (31)
Each subspace F n includes the previous ones, F n-1 ⊂ F n , and therefore inequalities

D 0 ≥ D 1 ≥ D 2 ≥ • • • . ( 32
)
The series D n converges to the actual diffusivity:

lim n→∞ D n = D. (33)
The subspace F 0 consists of constant functions. Therefore all sums in (30) vanish and the diffusivity in the zeroth order approximation reads This equation is valid for k-GEPs with any k; the same hold for all previous results of this section. In the rest of this work, we consider the simplest 2-GEP models if not stated otherwise.

D 0 = p τ0τ1 χ . (34)
For the 2-GEP, the zeroth order approximation (34) for the diffusion coefficient becomes

D 0 = p 10 + 2p 11 λ + ap 21 λ 2 1 + 4aλ + aλ 2 . ( 35
)
It is also possible to obtain an explicit form of D 0 as a function of ρ by substituting ( 22) into (35). We present D 0 in two cases. Specializing (35) to the particle-uniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF] we obtain

D 0 = 2(1 + λ) 2 2 + 4λ + λ 2 = 1 + ρ + 1 + 2ρ -ρ 2 2 1 + 2ρ -ρ 2 (36)
while in the site-uniform case [START_REF] Derrida | [END_REF] [23]

D 0 = (1 + λ) 2 1 + 4λ + λ 2 = 2 + 1 + 6ρ -3ρ 2 3 1 + 6ρ -3ρ 2 . ( 37
)
There is a symmetry between particles and vacancies for the 2-GEP with site-uniform rates. This symmetry implies the mirror symmetry of the diffusivity, that is, the invariance with respect to the ρ ↔ 2 -ρ transformation:

D(ρ) = D(2 -ρ) . ( 38
)
In terms of the fugacity, the mirror symmetry reads

D(λ) = D(1/λ) . (39)
Our iterative scheme gives approximations of the diffusivity agreeing with mirror symmetry. This is evident in the zeroth order, Eq. ( 37), and holds in all higher orders.

C. The first-order approximation

The case of n = 1 is more instructive, although we shall see that the final result is the same as for n = 0. When n = 1, the functions f depend only on τ 1 , that is, FIG. 3. Illustration of the procedure for n = 2, when the local function f depends on two adjacent sites and only f, f 1 and f 2 contribute to (30).

f (τ ) = f (τ 1 ) := f τ1 . There are three independent values, f 0 , f 1 and f 2 , so F 1 is a three-dimensional space.

Only j = 0 and j = 1, see Fig. 2, contribute to the sum j in Eq. (30). Writing τ 0 = r and τ 1 = s, we have

Q(f )(τ ) = =±1 p rs ( -f r-1 + f r -f s+1 + f s ) 2 , ( 40
)
where the summation over = ±1 combines the two terms on the right-hand side of Eq. ( 30) into a single sum. Averaging (40) with respect to the equilibrium measure we obtain

Q(f ) = r=1,2 s=0,1 =±1 p rs W r W s ( -f r-1 + f r -f s+1 + f s ) 2 .
(41) By a straightforward calculation, one finds

Q(f ) = 2 p rs + (f 0 -2f 1 + f 2 ) 2 4λ 2 p 11 Z 2 . ( 42
)
This is minimized when f 0 -2f 1 + f 2 = 0 and gives D 1 = D 0 . Thus the first iteration does not lead to any improvement.

D. The second-order approximation

Now f (τ ) = f (τ 0 , τ 1 ) =: f τ0,τ1 can take nine different values. In the functional Q(f ), the sites j contributing to the sum j in Eq. ( 30) are j ∈ {0, 1, 2}, see Fig. 3. We write τ -1 = q, τ 0 = r, τ 1 = s, τ 2 = u. The average of Q(f ) becomes

Q(f ) = q,r,s,u∈{0,1,2} =±1 p rs W q W r W s W u (43) ( -f q,r-1 + f q,r -f r-1,s+1 + f r,s -f s+1,u + f s,u ) 2 .
To determine the minimum, we must solve ∂ Q(f ) ∂f ξη = 0 (ξ, η ∈ {0, 1, 2}). There are five independent equations characterizing the solution space. Four of them are homogeneous

f 20 -2f 11 = 0, f 10 -f 00 -f 11 = 0, f 22 + f 00 -2f 11 = 0, f 21 + f 12 + f 10 -f 20 -2f 11 = 0. (44)
An additional inhomogeneous equation is Using ( 44)-( 45) we find the minimum of Q(f ) , from which

f 12 + f 01 -f 02 -f 11 = p 21 + p 10 -p 20 V (
D 2 = D 0 - 2(p 21 + p 10 -p 20 ) 2 λ 2 (a -1 + 4λ + λ 2 )V Z . ( 46
)
The second order approximation gives a better prediction for the diffusion coefficient:

D 2 < D 1 = D 0 .
The only exception is the misanthrope process (13) which is gradient and therefore the zeroth order approximation is already an exact answer. Specializing the general expression (46) to the particleuniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF] we obtain

D 2 = 2(7 + 21λ + 23λ 2 + 13λ 3 + 3λ 4 ) (2 + 4λ + λ 2 )(7 + 7λ + 3λ 2 ) . ( 47
)
Figure 4 (a) shows D 0 (ρ), D 2 (ρ) and simulation results obtained for the system with open boundaries (Sec. IV). We see that simulation data almost perfectly match analytical results for D 2 .

For the site-uniform rates (12) we get

D 2 = 4 + 13λ + 16λ 2 + 13λ 3 + 4λ 4 (1 + 4λ + λ 2 )(4 + 5λ + 4λ 2 ) . (48)
The second-order approximation obeys the mirror symmetry (39) as expected.

E. The third and higher order approximations

In the third order, j ∈ {0, 1, 2, 3} contribute to the sum j in Eq. ( 30), see Fig. 5. We denote by τ -2 = , τ -1 = q, τ 0 = r, τ 1 = s, τ 2 = u, τ 3 = v the relevant occupation numbers and recast Q(f ) into

Q(f ) = ,q,r,s,u,v∈{0,1,2} =±1 p rs W W q W r W s W u W v ( -f ,q,r-1 + f ,q,r -f q,r-1,s+1 + f q,r,s -f r-1,s+1,u + f r,s,u -f s+1,u,v + f s,u,v ) 2 .
(49) We must solve ∂ Q(f ) ∂f ξηζ = 0 (ξ, η, ζ ∈ {0, 1, 2}). There are 17 independent equations which are listed in Appendix C.

Lengthy but straightforward calculations give

D 3 = a -1 A(λ) (a -1 + 4λ + λ 2 )B(λ) . ( 50
)
where A(λ) and B(λ) are polynomials in λ of the 8th and 6th degree, respectively. We do not display long explicit formulas for A and B for the general rates p rs and limit ourselves with two examples. FIG. 6. Difference between nth and (n -1)st order approximations vs ρ (a) and n (b). The plots are for the 2-GEP with particle-uniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF]. The quality of approximation is the worst around half-filling; in the examples of (b), the quality is worst when λ = 2.81 corresponding to ρ = 1.38.

For the particle-uniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF], we have a = 1 2 and A(λ) = 2870 + 16194λ + 40422λ 2 + 59250λ 3 + 56100λ 4 + 35274λ 5 + 14410λ 6 + 3499λ 7 + 385λ 8 , B(λ) = 2870 + 10454λ + 17064λ 2 + 15772λ 3 + 8707λ 4 + 2729λ 5 + 385λ 6 .

For the site-uniform rates [START_REF] Derrida | [END_REF], we have a = 1 and

A(λ) = 225 + 1428λ + 4140λ 2 + 7396λ 3 + 8886λ 4 + 7396λ 5 + 4140λ 6 + 1428λ 7 + 225λ 8 , B(λ) = 225 + 978λ + 2079λ 2 + 2620λ 3 + 2079λ 4 + 978λ 5 + 225λ 6 .
In Fig. 4 (b), we plot the curves D 0 , D 2 and D 3 as well as the simulation result around ρ = 1. We find that the difference between D 2 and D 3 is almost the same magnitude of the statistical error of the simulation.

It is possible to determine analytical expressions for D n for any finite n, but calculations become cumbersome as n increases. For n = 4, for instance, even in a simple case of particle-uniform rates the analytical expression for D n has the form (50), with A(λ) and B(λ) being polynomials with integer coefficients of degree 30 and 28, respectively (the coefficients are huge, some are of the order of 10 24 ). Figure 6 shows D n-1 -D n for n ≤ 6 for the case of the particle-uniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF]. The difference D 5 -D 4 , for instance, is less than 10 -5 in the entire density region 0 < ρ < 2. We expect that D n-1 -D n decreases faster than algebraically.

In a previous work [34], we calculated the diffusion coefficient for k-GEPs using an approach that actually gave the zeroth-order approximation, Eq. ( 34). This prediction was in good agreement with Monte Carlo simulations, but further investigations [35,36] of the 2-GEP revealed a discrepancy between the prediction and the actual value. The above approach shows how to improve this approximation in a systematic manner.

IV. GEPS WITH OPEN BOUNDARIES

On a ring, the 2-GEP and some classes of the k-GEPs (including particle-uniform and site-uniform cases) satisfy detailed balance, and equilibrium steady states are described by a product measure [START_REF] Cocozza-Thivent | [END_REF]. On an open chain connected to reservoirs at different densities, however, these processes are non-equilibrium and their steady states are generally unknown. (The simple exclusion process k = 1 is an exceptional case where the steady state is known, see e.g. [START_REF] Blythe | [END_REF][START_REF] Derrida | [END_REF].) In the large system size limit, the deviation from equilibrium is small and one can use Fick's law to describe the hydrodynamic behavior [START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF]. In this section we study numerically more subtle properties such as finite-size corrections to the product measure. Specifically, we consider the 2-GEP on a finite chain with L-1 bulk sites connected to reservoirs with densities ρ 0 and ρ L , see Fig. 7. The couplings to the reservoirs are described by injection and extraction rates α, β, γ, δ, see Fig. 7. These rates are determined by the densities ρ 0 , ρ L and the bulk hopping rates, see Appendix D for details. We use simulations and scaling considerations to shed light on non-equilibrium steady states. We perform simulations for the 2-GEP with particle-uniform rates. In our simulations, we consider extreme boundary densities, ρ 0 = 2 and ρ L = 0. Further, we perform a time average over 10 7 ≤ t ≤ 10 9 and an ensemble average over 10 independent runs.

A. Density profiles in the open system

In the long-time limit, the density profile in the open system becomes stationary and the current is uniform. The knowledge of D(ρ) allows us to determine the density profile by solving the stationary diffusion equation d dx D(ρ) dρ dx = 0 and by imposing the boundary conditions matching the densities of reservoirs, ρ(0) = ρ 0 and ρ(1) = ρ L . (Here, x = i/L is the rescaled length.) This gives the implicit form of the density profile

ρ(x) ρ0 D(ρ )dρ = x ρ L ρ0 D(ρ )dρ , (51)
valid for 0 ≤ x ≤ 1. Since this form contains both ρ 0 and ρ L , the density profile ρ(x) depends, of course, on these boundary densities. Replacing the exact diffusivity D(ρ) by D n (ρ) (for n = 1, 2, 3, . . .) derived in the previous section, leads to increasingly accurate approximations ρ n (x) of the true profile ρ(x). As shown in Fig. 8 (a), the difference between ρ(x) and ρ 0 (x) is already too small to be visible when L = 1024, namely it is less than 0.004 in absolute value, see Fig. 8 (b). However, the discrepancy between ρ(x) and ρ 0 (x) is systematic and is not due to statistical errors: Fig. 8 (c) indicates that it does not vanish in the L → ∞ limit. Figure 8 (b,c) also demonstrates that ρ(x) -ρ 0 (x) is well fitted by ρ 3 (x) -ρ 0 (x). Generally, ρ n (x) with increasing n provide more and more accurate predictions for ρ(x).

If the gradient condition is satisfied, as for the misanthrope process, the exact bulk diffusivity is given by D 0 (ρ) and therefore ρ 0 (x) becomes identical to ρ(x) in the hydrodynamic limit L → ∞.

B. Finite-size corrections to the product measure

Because the stationary measure is not factorized in finite systems with open boundaries, the knowledge of the local densities

ρ i = τ i = P[τ i = 1] + 2P[τ i = 2] (52)
is not sufficient to determine the stationary state. In an open system of size L, deviations from the product measure scale as L -1 in the leading order. In order to probe these deviations, we first consider the quantity

Γ i = aP[τ i = 1] 2 -P[τ i = 0] P[τ i = 2] (53)
which vanishes for infinite systems and for finite systems on the ring, see Eq. (18). Figure 9 shows simulation results in the system with particle-uniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF] and reservoir densities ρ 0 = 2 and ρ L = 0. The plots indicate that Γ i acquires a scaling form:

Γ i L -1 ω(x) (54)
in the scaling limit The difference ρ 0 (x) -ρ(x) for x = 3/8 and 7/8, as a function of system size. The difference ρ 0 (x) -ρ 3 (x), for the same two values of x, are also shown. For the three panels, we choose the boundary densities (ρ0, ρL) = (2, 0).

L → ∞, i → ∞, x = i L = finite. ( 55
Next we consider the one-point functions P[τ i = r]. The corresponding finite-size corrections have the form similar to (54):

P[τ i = r] -W r (ρ i ) L -1 κ r (x) (56)
with scaling functions κ r obeying

κ 0 = - κ 1 2 = κ 2 = - ω(x) 1 -(1 -4a)ρ(x)[2 -ρ(x)] . (57)
We verified (56) for the particle-uniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF] and expect (56) to hold for a generic choice of hopping rates. The pair correlation function between two adjacent sites plays a crucial role in the following analysis. We use notations

X rs i = P[τ i = r ∧ τ i+1 = s] (58a) Y rs i = P[τ i = r] P[τ i+1 = s] (58b)
with r, s ∈ {0, 1, 2} and focus on the connected version of the pair correlation function

C rs i = X rs i -Y rs i . (59)
For a periodic or infinite system, C rs i = 0. However, this is generally not true for open systems of size L. In the scaling limit (55), the pair correlation function acquires a scaling form

C rs i L -1 ϕ rs (x) . ( 60
)
Note that these correlation functions satisfy simple sum rules ϕ r0 (x) + ϕ r1 (x) + ϕ r2 (x) = 0 (r = 0, 1, 2) (61a)

ϕ 0s (x) + ϕ 1s (x) + ϕ 2s (x) = 0 (s = 0, 1 , 2) (61b)
obtained by summing over the possible occupations of the right-hand side of (59). Only four of the nine functions ϕ rs are independent. In Fig. 10, all nine functions are plotted. The sum rules are obeyed thereby providing a check of our simulations.

C. Hydrodynamic limit for the current

We finally derive an alternative formula for the diffusion coefficient D(ρ) by taking the continuous limit of the exact expression for the microscopic current. It is crucial to keep the dominant corrections to the product measure defined in (60). This approach allows us to identify the missing contribution in the naive hydrodynamic limit of [34]. The alternative formula for D(ρ) is again checked by numerical simulations.

The current between site i and i + 1 reads

J i = 1≤r≤2 0≤s≤1 p rs (X rs i -X sr i ) (62)
and it can be re-written as

J i = J i + (p 21 + p 10 -p 20 )(C 21 i -C 12 i ) , (63a)
J i = 1≤r≤2 0≤s≤1 p rs (Y rs i -Y sr i ), (63b)
see Appendix E for details. The term J i involves only the products Y rs i of one-point functions. The second term on the right-hand side of (63a) vanishes when the 2-GEP satisfies the gradient condition (i.e., for the misanthrope process).

In the scaling limit (55), the term J i becomes

J(x) = - 1 L D 0 (ρ) dρ dx , (64)
where D 0 is the zeroth order approximation of the diffusivity (see Appendix F for details). Note that (64) is identical to the naive hydrodynamic limit [34].

The second term on the right-hand side of (63a) also simplifies in the scaling limit (55), namely it turns into

L -1 µ(x) with µ(x) = (p 21 + p 10 -p 20 ) ϕ 21 (x) -ϕ 12 (x) . (65)
Combining ( 64) and (65) we arrive at

LJ(x) -D 0 (ρ) dρ dx + µ(x), (66)
which in conjunction with Fick's law implies that

D(ρ) = D 0 (ρ) -µ(x) dx dρ . (67)
Here x = x(ρ) is the inverse function of the stationary density profile ρ(x) given by (51). Equation (67) allows us to obtain numerical plots of D(ρ) since x = x(ρ) and µ(x) can be determined accurately from simulations.

(This procedure was used in preparing Fig. 4.) Rewriting Eq. (67) and using the implicit equation for the profile (51) we obtain

µ(x) = D 0 (ρ) -D(ρ) dρ dx = 1 - D 0 (ρ) D(ρ) I (68)
with

I ≡ - ρ L ρ0 D(ρ)dρ. (69)
Replacing D by D n , one obtains Simulation results for µ(x) are well fitted by µ 3 (x), see Fig. 11. In Fig. 11 we also plot I defined in (69) as a function of L. We observe an excellent agreement between I (simulations) and I 3 (theory) as L increases. We also verify that I ≡ -ρ L ρ0 D 0 (ρ)dρ , converges to I 0 when L → ∞. This confirms the usefulness of Eq. (68) for numerical measurements of the diffusion coefficient, as an alternative to the standard Green-Kubo formula.

µ n (x) = 1 - D 0 (ρ) D n (ρ) I n , I n = - ρ L ρ0 D n (ρ)dρ . ( 70
V. CONCLUSIONS Interacting many-particle systems generically exhibit macroscopic behaviors that can be described by fluctuating hydrodynamics. This is particularly well understood in the realm of stochastic lattice gases with symmetric hopping. In these diffusive lattice gases, the hydrodynamic behavior is described by the single scalar function, the coarse-grained density ρ(x, t), satisfying the diffusion equation. Fluctuations and large deviations in diffusive lattice gases are described by two coupled scalar fields evolving according to equations of the macroscopic fluctuation theory. Remarkably, all details of the microscopic hopping rules are encapsulated in the diffusion coefficient D(ρ) and the conductivity σ(ρ). These transport coefficients enter into the diffusion equation and the governing equations of the macroscopic fluctuation theory. The determination of the transport coefficients is a challenge even for the simplest diffusive lattice gases. Apart from very special models, known as gradient systems, for which a simple perturbation approach provides the correct answer, there exists no general technique to derive closed formulas for the transport coefficients.

In this work, we used the Varadhan-Spohn variational formula for the diffusion coefficient and employed an ap-proach resembling the Ritz method to derive analytical approximations for D to an arbitrary degree of precision. The Varadhan-Spohn formula is essentially the Green-Kubo formula for diffusive lattice gases. We illustrated our approach by investigating the 2-GEP, a generalization of the symmetric exclusion process in which a site can accommodate at most two particles. The 2-GEP is nongradient gas with unknown transport coefficients. For the 2-GEP in one dimension we showed that the Varadhan-Spohn variational formula can be used to derive analytical approximations to D in a systematic manner-after a few iterations, the precision of the order of a few parts per million is reached. Thus the simplest approximations are already remarkably accurate resembling the success of continuous fraction approximations of irrational numbers. We emphasize that our approximate expressions for the diffusion coefficient are intrinsically different from perturbative, Taylor-type, expansions around some special value of the density; on the contrary, they are global and uniform estimates, valid on the whole range of admissible densities. In order to check for the practicality of our approximation procedure, we performed simulations in open systems coupled to reservoirs at different densities. The precise variational estimate of the diffusivity allowed us to study the density profile in the open system. In addition, an analysis of the finite-size corrections to the equilibrium measure led to an alternative formula for the diffusivity that can be used for high precision numerical measurements.

Nongradient lattice gases are the rule, not the exception, and the procedure used in this work can be extended in various directions. One could study one-dimensional processes with local hopping rates that depend in an arbitrary way on the occupation numbers of the original and of the target site. Further, the variational formula is valid in arbitrary spatial dimension. For instance, one can apply it to study the diffusion coefficient and conductivity of kinetically constrained lattice gases [START_REF] Sollich | Reviews of kinetically-constrained lattice gases: F. Ritort[END_REF] in two and higher dimensions. Different variational methods have been applied to a number of lattice gas models in Refs. [START_REF] Za Luska-Kotur | [END_REF], while techniques similar to ours have been recently introduced to approximate the thermal conductivity for stochastic energy exchange models [39,40]. It would be interesting to clarify connections between these approaches. Another avenue for further research concerns the phenomenon of self-diffusion describing the diffusivity of a tagged particle in a lattice gas. The coefficient of self-diffusion is unknown even for diffusive lattice gases satisfying the gradient condition, e.g. for the simple exclusion process in two and higher dimensions. The coefficient of selfdiffusion can be written in a variational form [START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF], and it is perhaps possible to derive excellent analytical approximations for this coefficient using a procedure similar to the one developed in this article. density of defects (0's and 2's) is known to decay as t -1/4 . (This is certainly valid when both species of defects, 0's and 2's, have the same diffusion coefficients, that is when p 10 = p 21 .) Thus at the half-filling, ρ = 1, the 2-GEP with p 11 = 0 algebraically approaches to the uniform final state where each site is occupied by one particle.

FIG. 2 .

 2 FIG.2. Illustration of the procedure of computing Q(f ) for n = 1 when the local function f (τ ) = f (τ1) depends only on a single site and hence only f 0 ≡ f and f 1 contribute to(30).

  45) with V = p 21 + p 10 + p 20 + p 21 + 2p 11 λ + p 10 aλ 2 Z .

FIG. 4 .FIG. 5 .

 45 FIG. 4. (a) Zeroth-and second-order approximations of the diffusivity over the entire density range. (b) Zeroth-, second, and third-order approximations in the density range 0.9 ≤ ρ ≤ 1.1. Numerical results (dotted lines) were obtained by simulating an open system of length L = 1024.

1 FIG. 7 .

 17 FIG. 7.Illustration of the 2-GEP with open boundaries. The "sites" 0 and L correspond to reservoirs with (in general) different densities. These boundary conditions are realized by transitions of the states at sites 1 and L -1, which are indicated by arrows.

FIG. 8 .

 8 FIG. 8. (a) ρ(x) for the system of size L = 1024 and ρ 0 (x) obtained from Eq. (51) with replacement D → D0. (b) The difference between ρ 0 (x) and ρ(x) observed in simulations. The difference ρ 0 (x) -ρ 3 (x) is shown for comparison. (c)The difference ρ 0 (x) -ρ(x) for x = 3/8 and 7/8, as a function of system size. The difference ρ 0 (x) -ρ 3 (x), for the same two values of x, are also shown. For the three panels, we choose the boundary densities (ρ0, ρL) = (2, 0).

FIG. 9 .

 9 FIG.9. The scaled correlation function LΓi vs the scaled coordinate i/L for three different system sizes. The data collapse supports the emergence of the scaling expressed by (54).

FIG. 10 .

 10 FIG. 10. Correlation functions (59) multiplied by the system size L. The data collapse supports the emergence of the scaling behavior (60). The sum rules (61a)-(61b) are obeyed.

FIG. 11 2 0 2 0

 1122 FIG. 11. (a) L(C 21 i -C 12 i ) vs x = i/L for three different system sizes; µ3 [Eq. (70)] is shown for comparison. (b) Simulation results for I = 2 0 D(ρ)dρ and I = 2 0 D0(ρ)dρ . The estimates In with n = 0, 2, 3 are also shown.

ACKNOWLEDGMENTS

CA thanks Isfahan University of Technology, Institute for Research in Fundamental Sciences (IPM), and IPhT, CEA Saclay for hospitality. The work of PLK was partially supported by the BSF Grant No. 2012145. PLK is also grateful to IPhT, CEA Saclay for excellent working conditions, and to Grant Number ANR-10-LABX-0039-PALM for supporting the visits. We are thankful to S. Mallick for a careful reading of the manuscript.

Appendix A: Extreme versions of the 2-GEP For the 2-GEP, the expressions (14) for the one-site probabilities indicate that unusual behaviors may occur in the extreme cases of a = 0 and a = ∞. The latter case, a = ∞, is realized when p 20 = 0. The 2-GEP in this case is not a diffusive lattice gas. Instead, starting from an arbitrary initial configuration, the system gets trapped in a jammed configuration like . . . 002020002222022202022200002022202 . . . where each site is empty or occupied by two particles. In the long-time limit the configuration looks like . . . 020201022221222020222001020222102 . . . These 1's diffuse, 10 ↔ 01 and 12 ↔ 21, and annihilate, 11 → 20 or 11 → 02, i.e., 1's disappear due to the singlespecies diffusion-controlled annihilation (see [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF] for a review). In one dimension, the density of 1's decays as t -1/2 . Thus, the 2-GEP with p 20 = 0 is not diffusive and the system gets trapped in a jammed final state.

The case of a = 0 is realized when p 11 = 0 (see also [29]). Note that, in this case, D(ρ) is not continuous at ρ = 1, and correspondingly d 2 F dρ 2 diverges at this point. Elsewhere the Einstein relation (2) should be satisfied. If the global density is ρ < 1 and additionally we start with a configuration where each site is occupied by at most one particle, τ i (0) ≤ 1, such occupation arrangement will persist, so we recover the classic exclusion process. Similarly if the global density is ρ > 1 and additionally we start with a configuration where each site is occupied by at least one particle, τ i (0) ≥ 1, the process can be mapped into the classic exclusion process after interpreting a site with two particles as occupied and a site with one particle as empty. Generally when ρ = 1, the 2-GEP with p 11 = 0 essentially reduces to the classic exclusion process after an earlier regime when the system reaches state with τ i ≤ 1 (if ρ < 1) or τ i ≥ 1 (if ρ > 1). A novel behavior can only occur when ρ = 1. In the long time, the configuration will look like . . . 11111011111111112111110111111111211111 . . . with 0's and 2's diffusing and annihilating upon colliding: 02 → 11, 20 → 11. For this two-species diffusioncontrolled annihilation process (see [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF] for a review), the Appendix B: Validity of the product measure for the 3-GEP

Here we briefly consider the 3-GEP to appreciate possible qualitative differences with the 2-GEP. The 3-GEP is generally characterized by nine non-vanishing rates p rs with r = 1, 2, 3 and s = 0, 1, 2. The detailed-balance relations ( 17) lead again to (18) and additionally to

The equilibrium probabilities take the form (5), viz.

with a 3 = a 2 p 12 /p 30 and the same a 2 = p 11 /p 20 as for the 2-GEP. All this is valid, however, only when the constraint

following from (B1) is obeyed. Thus for the 3-GEP, the factorization holds when the constraint (B2) is satisfied; similar conditions for the k-GEPs with k ≥ 4 are discussed in Ref. [START_REF] Cocozza-Thivent | [END_REF].

Appendix C: Details for n = 3

Solving ∂Q ∂f ξηζ = 0 (ξ, η, ζ ∈ {0, 1, 2}), one arrives at 13 homogeneous relations (we use the shorthand notation

and 4 inhomogeneous relations

Here U i = A i Z/B, A i are polynomials of the 4th degree in λ and B is the polynomial of the 6th degree. For the particle-uniform rates [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF],

For the site-uniform rates [START_REF] Derrida | [END_REF],

For the misanthrope process A i = 0. The expressions for B(λ) for the particle-uniform and site-uniform rates are presented in subsection III E.

Appendix D: Microscopic coupling to the reservoirs At the left boundary, the site number 1 of the open system is connected to a reservoir at density ρ 0 . Particles are injected to site 1 at rate α(τ 1 ) and removed from site 1 at rate γ(τ 1 ). These rates depend on the occupation τ 1 of site 1, so there are four different values: α(0), α(1), γ(1) and γ(2). Imposing a local equilibrium condition with the left reservoir at density ρ 0 implies the following constraints:

The weights W 0 , W 1 and W 2 are the 2-GEP equilibrium probabilities ( 14) at density ρ 0 , the corresponding fugacity λ 0 = λ(ρ 0 ) is obtained from (22). Relations (D1) leave some freedom in choosing the boundary rates. The following choice is suitable:

The local detailed balance relations are readily checked:

At the right boundary, the system is at local equilibrium with a reservoir at density ρ L ; particles are extracted from site L -1 at rate β(τ L-1 ) and injected with rate δ(τ L-1 ). A suitable choice of boundary rates is

where now the weights W s are evaluated for the value of λ L , obtained from ρ L using (22).

In general, the product measure does not give the correct stationary state, except for the case ρ 0 = ρ L . The difference between the reservoir densities ρ 0 and ρ L gives a finite current even in the stationary state.

Simulations shown in this work were performed in the situation when the boundary densities are extreme, viz. the left reservoir is fully packed, ρ 0 = 2, while the right reservoir is empty, ρ L = 0; the corresponding values of the fugacity are λ = ∞ and λ = 0, respectively. Accordingly the injection and extraction rates are

First, we show that Eq. ( 62) is equivalent to

Note that the top line on the right-hand side of Eq. (E1) is in a gradient form, but the bottom line is not (unless p 10 -p 20 + p 21 = 0). We use the following identities

for τ i ∈ {0, 1, 2}. The Kronecker delta symbol δ τi,r indicates if site i is occupied by r particles. Similar expressions are readily written for two-point probabilities

The expression (62) for the current J i is thus given by the expectation value

Introducing the observables

we obtain

leading to (E1). Performing analogous calculations for J i defined in (63b), we find that J i -J i is given by

times (p 21 +p 10 -p 20 ). One can verify that (E2) is identical to C 21 i -C 12 i . This completes the proof of (63a)-(63b).

Appendix F: Continuous limit of Ji

We now derive (64), the continuous limit of J i . We use the scaling form (56) and obtain We emphasize that the finite-size corrections κ r cancel out in (F3). Substituting ( 14) into (F3) we obtain

p rs W r W s which indeed reduces to (64) once we recall that D 0 is given by (34).


Dan Eude Kpannieu

Martine Mallet

Mustapha Lacina Coulibaly

Mustapha Abdelmoula

Christian Ruby
email: christian.ruby@univ-lorraine.fr

Phosphate removal from water by naturally occurring shale

Keywords: phosphate, adsorption, iron, geomaterials, Ivory Coast

come

Introduction

Phosphorous is an essential nutrient for the growth of aquatic plants and algae but its excessive supply in water bodies causes eutrophication resulting in increased growth of algae, the depletion of oxygen and the disruption in the balance of organisms. Excessive phosphate inputs to waters originates from anthropogenic activities such as wastewaters discharges and agricultural run-off. In wastewater treatment technology, chemical and biological removal techniques are the most effective and well-established methods but their high cost and complexity of implementation prevent their application both in rural areas and in developing countries. Recently, attention has been focused on developing alternative methods for phosphate removal from water. Among these methods, adsorption appears as an attractive option, owing to its low cost and easy application [START_REF] Boyer | Comparison of low-cost and engineered materials for phosphorus removal from organic-rich surface water[END_REF][START_REF] Hongbin | Simultaneous removal of ammonium and phosphate from eutrophic waters using natural calcium-rich attapulgite-based versatile adsorbent[END_REF].

A large variety of materials with the availability to remove phosphate were thus described in the literature [START_REF] Johansson-Westholm | Review-Substrates for phosphorus removal-Potential benefits for on-site wastewater treatment?[END_REF][START_REF] Vohla | Filter materials for phosphorus removal from wastewater in treatment wetlands-A review[END_REF][START_REF] Wendling | Nutrient and dissolved organic carbon removal from natural waters using industrial by-products[END_REF]. High sorption capacity, low cost, availability in large quantity and environmentally-friendly nature are among the most important parameters in the choice of the sorbent. The present paper aimed to develop a material satisfying all these criteria with the final objective of applying this optimized material in the treatment of wastewater released into the Ebrié lagoon in Ivory Coast.

In fact, the eutrophication of the Ebrié lagoon is dramatic in large areas, making the lagoon actually unsafe for habitation sites, for the practice of lucrative activities such as bathing and fishing [START_REF] Briton Bi | Evaluation of the Abidjan lagoon pollution[END_REF][START_REF] Koffi | Pollution of a tropical lagoon by the determination of organochlorine coumpounds[END_REF][START_REF] Tuo | Assessment of waters contamination by nutrients and heavy metals in the Ebrie Lagoon (Abidjan, Ivory Coast)[END_REF][START_REF] Yao | La pollution des eaux de la zone urbaine d'une lagune tropicale par les matières oxydables (lagune Ebrié, Côte d'Ivoire)[END_REF]. In this context, shale, sandstone and laterite are natural geological materials easily available in large quantities in Ivory Coast. These low cost and environmentally friendly materials contain iron, calcium and aluminum and are thus expected to be promising phosphate adsorbents.

The main objective of the present study was to assess and compare the removal efficiencies of these three materials. Batch adsorption experiments were carried out as a function of adsorbent dosage, pH and contact time. The adsorption kinetics and Langmuir and Freundlich isotherms were also studied. Then, the performance of the three materials for phosphate removal were evaluated in flow through column experiments. The originality of the present work lied in the understanding of the relation existing between the physicochemical properties of the adsorbents and their ability to remove phosphate from water. A particular attention was devoted to quantify the dissolution properties of the minerals (release of Fe II , Al III ,, Mg II , Ca II in solution) in order to estimate the respective roles of dissolution precipitation and adsorption processes in the phosphate removal from water. This study is a necessary preliminary work before using the best determined adsorbents in a constructed wetland system.

Materials and methods

Chemicals

All chemicals used were analytical grade. All experiments were conducted in polypropylene vessels that were soaked in 36% hydrochloric acid (HCl) before use and rinsed with deionized water (Milli-Q, 18.2 M cm -1 ) and ethanol. Deionised water was used for preparing all solutions.

Adsorbent.

Shale, sandstone and laterite were collected in three regions of Ivory Coast, i.e. shale from the Center region (Toumodi-Lomo North), sandstone from the south region (Abidjan-Akouedo) and laterite from the Eastern region (Abengourou-Sankadjo). Samples were prepared using the same following procedure. Large blocks extracted from the soil were first coarsely screened, washed with deionised water to remove the surface-adhered particles and dried in an oven at 50°C for 24 hours. The samples were then crushed again and sieved (AFNOR sieves) to obtain a granular size  400 μm. Powder samples were finally washed again and dried at 70°C for 24 hours before their use. The chemical composition of the three substrates (Table 1) revealed that all three materials were primarily a mixture of oxides, the major one being SiO2, Al2O3 and Fe2O3. The relative amount of these oxides differed from one material to the other. Sandstone was iron-rich while shale and laterite were Al-rich materials. All three subtrates contained large amount of Si. The specific surface area of the three powders materials (  400 µm) was determined by multipoint N2-BET analysis using a Belsorp max surface area analyser (Bel, Japan). The specific surface area of shale, sandstone and laterite were 6.3, 28.5 and 37.6 m 2 g - 1 respectively.

Adsorption studies.

Phosphate adsorption onto the three substrates was studied in both batch and column experiments at room temperature. A stock phosphate solution (250 mg L -1 ) of phosphate was prepared by using sodium dihydrogen phosphate (H2NaO4P.2H2O) and diluted in subsequent experiments.

Batch experiments.

The effect of sorbent dose was examined in a series of experiments that used the same initial phosphate concentration (25 mg PO4 L -1 ) while using different adsorbent doses between 0 and 160 g L -1 . 50 mL flasks were used. The pH of the suspension was maintained at a defined value of 6.5-7.0 by manually adding 0.5 M HCl or 0.5 M NaOH at regular intervals. Flasks were capped and stirred magnetically at 300 rpm for 24 hours to ensure approximate equilibrium. At the end of this adsorption period 10 mL were withdrawn from the reaction flasks and filtered through 0.22 µm polypropylene syringe filters for phosphate and aqueous metal species (Ca, Mg, Al and Fe) analyses by ICP-AES. All experiments were carried out in duplicate to check the reproducibility of the experiments. The pH effect on phosphate adsorption was examined with a similar procedure using a 25 mg L -1 initial phosphate concentration and 80 g L -1 (shale) and 80 and 8 g L -1 (sandstone and laterite) adsorbent dosages while maintaining pH at different values in the range 2-12. Phosphate sorption kinetics were examined in 1 L Erlenmeyer flasks, using an initial phosphate concentration of 25 mg L -1 , various adsorbent doses, and for ~ 60 h.

During this reaction period, samples were withdrawn from the flasks at regular intervals for phosphate and aqueous metal species determination. The pH of the suspension was maintained at a defined value (6.5-7.0). Sorption isotherms experiments were carried out at pH = 6.5-7.0 in 50 mL Erlenmeyer flasks with varying initial concentrations of phosphate ions (5-250 mg L - 1 ) and a constant adsorbent dosage of 80 g L -1 (shale) and 8 g L -1 (sandstone and laterite). The removal efficiency R (%) of phosphate onto adsorbents and phosphate uptake qt (mg g -1 ) at time

t were calculated from Eqs 1-2: i C ) t C i C ( R   (1) m V ) t C i C ( t q    (2)
where Ci and Ct are the initial and at time t phosphate concentration (mg L -1 ) respectively; V is the volume of the solution (L) and m is the mass of material (g). Note that phosphate uptake at near equilibrium conditions (qe, mg g -1 ) is obtained at the equilibrium phosphate concentration Ct = Ce.

Column experiments.

Plexiglas columns of 500 mm length and 37 mm internal diameter corresponding to a volume Vtot of 473 mL were used. The columns were packed with dried shale, sandstone or laterite (grain size of 1-2 mm) corresponding to a total mass of 593, 661 and 679 g, respectively. The porosity of all substrates was very close from 50 %; indeed the porous volume VP for the three materials was very similar, i.e. VP = 236, 246 and 241 mL for shale, sandstone and laterite,

respectively. An upwards water flow direction was used with the help of a peristaltic pump (Watson Marlow) ensuring by this way water-saturated conditions. All the columns were first pre-equilibrated with demineralized water during 24 hours before performing the phosphate removal experiments. The pH of the phosphate solution ([PO4] = 25 mg L -1 ) introduced into the column was close to 5 and the pH measured in the outflow was relatively constant independently of the nature of the materials used in the columns with a mean pH value of 6.1, 6.0 and 6.4 for shale, sandstone and laterite, respectively. For both sandstone and laterite columns, a continuous flow of 1 mL min -1 with a corresponding Hydraulic Retention Time (HRT) of  4 hours was used. The total duration of the experiments was  2.5 months. In order to mimic situations that can be encountered in constructed wetlands where water supply is often provided discontinuously, a water flow circulating discontinuously was used for the shale column as described in our previous work [START_REF] Kpannieu | Removal of phosphate by shale from Ivory Coast in homogeneous reactor and hydrodynamic conditions: influence of soluble species[END_REF]. The water flow was 0.23 mL min -1 during a first period of 104 hours with a corresponding HRT  17 hours followed by a shutdown of the water pump corresponding to a HRT of  64 hours. The sequence (104 hours, 64 hours) was repeated 9 times for a total duration of the experiments of  2 months. At appropriate time intervals, 10 mL of solution was sampled in the column outflow and filtered through 0.22 µm polypropylene syringe filters before phosphate and metal species were analyzed with ICP-AES.

Characterization methods

The structural properties of shale, sandstone and laterite were investigated with X-Ray Diffraction (XRD) and Mössbauer spectrometry. XRD patterns were acquired with a PANalytical X'Pert Pro Multipurpose diffractometer using Co K radiation ( = 1.78026 Å; 40 kV, 40 mA). Diffraction patterns were obtained over a range of 10-85 o 2θ range with a step size of 0.036° and nominal time per step of 3 s. Transmission Mössbauer spectrometry was performed with a 50 mCi source of 57 Co in rhodium to identify and quantify iron phases. The Mössbauer spectra were calibrated relatively to the reference spectrum of a 25 µm thick pure iron foil recorded at room temperature. A cryostat equipped with a closed cycle helium station allows variable temperature measurement between 8 and 300 K. The spectra were fitted using

Lorentzian shape lines and, in some cases, Voigt profile analyses, a convolution of a Gaussian distribution with a Lorentzian shape line, was used. Residual phosphate (PO4-P) and metal species (Ca, Mg, Al and Fe) concentrations in solution were determined by Inductively Coupled Plasma-Atomic Emission Spectrometry (ICP-AES) using an Ultima spectrometer (HORIBA Jobin Yvon).

Results and discussion

3.1. Physicochemical properties of the adsorbents.

X-Ray Diffraction (XRD).

All three diffraction patterns (Fig. 1) displayed strong reflection of quartz (SiO2). Additional reflection lines of silicates minerals such as muscovite (KAl2(Si, Al)4O10(OH)2, albite (Na, Ca) (Si, Al)4O8) and clinochlore (Mg, Fe)5(Si3Al)O10(OH)8 were observed for shale while iron oxides such as goethite (-FeOOH) and hematite (Fe2O3) were observed for sandstone and laterite. Such mineral phases agreed relatively well with chemical composition reported in Table 1 except the absence of iron bearing minerals in shale. Because iron is well known to play a key role in phosphate adsorption [START_REF] Altundoğan | Removal of phosphates from aqueous solutions by using bauxite. I: Effect of pH on the adsorption of various phosphates[END_REF][START_REF] Mallet | Investigation of phosphate adsorption onto ferrihydrite by X-ray Photoelectron Spectroscopy[END_REF] a closer examination of iron phase minerals contained in the three substrates was carried out using Mössbauer spectrometry.

Mössbauer spectrometry.

The Mössbauer spectra of sandstone and laterite recorded at room temperature (Fig. 2) exhibited both two magnetic sextets attributed to hematite-Fe2O3 and goethite -FeOOH and a ferric paramagnetic doublet D1 (Table 2), in agreement with XRD analyses. Mössbauer analysis at lower temperature was needed to distinguish whether the paramagnetic doublet D1 recorded at room temperature corresponded to iron oxides or to Fe III species embedded into the structure of clays. The Mössbauer spectra of sandstone and laterite recorded at 9 K clearly showed that the doublet observed at 300 K has been transformed into a magnetic sextet attributed to -FeOOH (Fig. 3) [START_REF] Kone | Reductive transformation and mineralization of an azo dye by hydroxysulphate green rust preceding oxidation using H2O2 at neutral pH[END_REF]. The corresponding Mössbauer hyperfine parameters (Table 3) revealed that only a very low intensity (RA of  3 %) paramagnetic doublet remained on the 9 K Mössbauer spectrum of laterite and could be attributed to structural Fe III species into clays.

Therefore, for both sandstone and laterite, Fe III species were essentially included into ferric oxides (-Fe2O3 and -FeOOH) and the ferric doublet recorded at room temperature (Fig. 2) was attributed to superparamagnetic goethite [START_REF] Janot | Caracterisation de kaolinites ferrifères par spectrométrie Mössbauer[END_REF]. This meant that the size of some of the goethite particles in sandstone and laterite was quite small. This must be particularly true for laterite since the ferric doublet recorded at room temperature was relatively intense (RA of  55 %) and narrow and did not exhibit any relaxation effect (Table 2). This last result may be correlated to the BET analysis showing that laterite had the highest surface specific area.

The status of iron in shale was very different and the Mössbauer spectrum recorded at 12 K (Fig. 3) exhibited a ferric (D1) and a ferrous (D2) paramagnetic doublet. The total relative area of both doublets are relatively high representing  2/3 of the total spectra. This meant that for the shale most of the iron species were embedded into the clay minerals on the contrary to sandstone and laterite where iron was present into iron oxides. However, the shale did also contain some iron in goethite as shown by the magnetic sextet of -FeOOH observed at 12 K (Fig. 3).

Phosphate sorption in batch experiments

3.2.1. Effect of adsorbent dosage.

All three substrates displayed not surprisingly an increase in the removal efficiency with increasing adsorbent dosage (Fig. 4). A high removal efficiency of 100% was achieved with laterite and sandstone at 20 and 40 g L -1 respectively while removal efficiency remained less than 100 % in the whole sorbent dose studied for shale (Fig. 4a). In addition, the phosphate uptake strongly decreased as the dosage increased and reached a saturation level at high sorbent doses. Such results reflected that adsorbed phosphate species were distributed among more available sorption sites at high sorbent doses, thus decreasing the sorption capacity [START_REF] Ádám | Phosphorus retention in the filter materials shellsand and filtralite P® -Batch and column experiment with synthetic P solution and secondary wastewater[END_REF][START_REF] Babatunde | Characterization of aluminium-based water treatment residual for potential phosphorus removal in engineered wetlands[END_REF].

Effect of pH on phosphate adsorption.

The pH is an important factor that strongly influences the phosphate adsorption mechanisms at the solid-solution interfaces.. Small amounts of phosphate are determined at pH values lower than 8 and 10 for sandstone and laterite respectively which indicated a phosphate removal efficiency close to 100 % (Fig. 5). Beyond these two pH values, phosphate concentration strongly increased. A different behaviour was observed for shale. In fact, an increase in phosphate concentration was observed in the pH range of 2-11 followed by a strong decrease beyond pH 11. The insert in Figure 5 displays the phosphate concentration determined in solution for a 8 g L -1 laterite and sandstone dosage. The phosphate concentration increased onto the whole pH range studied. It can thus be concluded that increasing the pH resulted in a decrease of phosphate removal efficiency (at phosphate removal efficiencies lower than 100 %) except for shale for which a stronger removal efficiency was observed at pH 12. With the exception of this high pH value, the phosphate removal efficiency of shale was much smaller than those for sandstone and laterite.

To further investigate the phosphate removal mechanism, aqueous ions release was studied as a function of pH (Fig. 6). Strong release of Ca, Mg and Al, and to a lesser extent, Fe was observed for shale in the low pH range and in comparison to sandstone and laterite. Dissolution processes are thus much more pronounced for shale. Aqueous ions concentrations were compared both in blank experiments and after adsorbents reaction with a 25 mg L -1 phosphate solution (Fig. 7). Interestingly sandstone and laterite clearly exhibited a different behaviour from shale. In fact, low amounts of Ca, Mg, Al and Fe were determined in solution after shale reaction with phosphate ions and whatever the pH was. Additionally, calcium and magnesium concentrations did not significantly differ in blank experiments and after phosphate reaction with sandstone and laterite even at a lower dosage of 8 g L -1 (data not shown).Some differences in ions concentrations were observed at low pH where the presence of phosphate ions induced an increase in Fe and Al dissolution processes in laterite and sandstone respectively. From Fig.

5 to 7, it clearly appeared that shale displayed a different sorption behaviour from those of sandstone and laterite. The lower influence of dissolution processes towards phosphate removal in case of sandstone and laterite combined to a much more efficient phosphate removal clearly suggested that phosphate sorption was dominated by dissolution/reprecipitation processes for shale and adsorption processes for sandstone and laterite.

Phosphate adsorption kinetics.

Adsorption kinetics were carried out at 80 g L -1 (shale) and 8 g L -1 (laterite and sandstone) sorbent doses (Fig. 8). Results obtained at additional sorbent doses for laterite and sandstone are displayed in supporting information (Figs S1 and S2). All curves displayed almost the same shape. The adsorption rate was rapid at the beginning of the experiments and then decreased until a pseudo plateau was reached. The first step may be due to phosphate adsorption at the external surface of the sorbents while the second step probably reflected some slow sorption process due to phosphate diffusion into the pores of the sorbent [START_REF] Sanyal | Chemistry of phosphorus transformation in soils[END_REF][START_REF] Yang | Dewatered alum sludge: a potential adsorbent for phosphorus removal[END_REF]. The strength of the adsorbent decreased in the following order: laterite  sandstone  shale. In particular, the observed maximal phosphate uptake determined from Figure 8 were 3.2, 2.8, and 0.23 mg g -1 for laterite, sandstone and shale respectively. The uptake capacities of laterite and sandstone were higher than those reported for other natural materials such as bauxite (0.6 mg g -1 ) [START_REF] Drizo | Physicochemical screening of phosphateremoving substrates for use in constructed wetland systems[END_REF], limestone (0.55 mg g -1 ) [START_REF] Drizo | Physicochemical screening of phosphateremoving substrates for use in constructed wetland systems[END_REF], zeolite (0.25 mg g -1 ) [START_REF] Drizo | Physicochemical screening of phosphateremoving substrates for use in constructed wetland systems[END_REF], dolomite (0.3 mg g -1 ) [START_REF] Pant | Phosphorus retention capacity of root bed media of subsurface flow constructed wetlands[END_REF] and apatite (0.41 mg g -1 ) (Bellier et al. 2006). Additionally, all these literature data were obtained at higher initial phosphate concentrations (40-150 mg L -1 ). The sorption capacity of shale was thus smaller than those of laterite and sandstone by more than a factor of 10. The common pseudofirst-order and pseudo-second order equations were used to modelize the experimental data using a linear fitting (eqs 3-4):

Pseudo-first-order equation:

303 . 2 t 1 k e q log ) t q e q ( log    (3) Pseudo-second-order equation: e q t 2 e q 2 k 1 t q t   ( 4
)
where qe (mg g -1 ) is the equilibrium adsorption amount, qt (mg g -1 ) is the amount of phosphate adsorbed at any time t and k1 (min -1 ) and k2 (g mg -1 min -1 ) are the pseudo first-order and pseudo second-order rate constants respectively.

The linear plots of phosphate sorption kinetics obtained and the corresponding kinetics parameters are given in Figure 9, Figure S3 and Table 4 respectively. The kinetic data could only be well described using the pseudo-second-order kinetic model. High correlation coefficients were in fact obtained and the equilibrium sorption amounts (qcal) closely match the experimental data while the pseudo first-order-model did not fit the data. Similar kinetics results

were reported in the literature for various phosphate adsorbents [START_REF] Cai | Competitive adsorption characteristics of fluoride and phosphate on calcinated Mg-Al-CO3 layered double hydroxides[END_REF][START_REF] Namasivayam | Recycling Fe(III)/Cr(III) hydroxide, an industrial solid waste for the removal of phosphate from water[END_REF][START_REF] Zhu | Adsorption of phosphate on hydroxyaluminum-and hydroxyiron-montmorillonite complexes[END_REF]. Such results indicated that the adsorption process of phosphate on these three materials may be based on chemisorption. As the sorbent dose increased and the other experimental conditions remained unchanged, the value of k2 increased.

3.2.5. Phosphate adsorption isotherms.

Isotherms are commonly used to describe the sorption behaviour of a given adsorbent [START_REF] Yang | Adsorption of phosphate from water on lake sediments amended with zirconium-modified zeolites in batch mode[END_REF][START_REF] Liu | Removal of phosphate anions using the modified chitosan beads: Adsorption kinetic, isotherm and mechanism studies[END_REF].. The sorption behaviour of laterite and sandstone was clearly different from that of shale (Fig. 10a). As it is the case in most studies, the shale curve displayed a saturation point at which increasing the initial phosphate concentration did not cause a significant increase in the phosphate sorption capacity [START_REF] Zhu | Adsorption of phosphate on hydroxyaluminum-and hydroxyiron-montmorillonite complexes[END_REF]. The presence of this saturation level allowed us fitting the data by using the common Langmuir and Freundlich models (Fig. 10b). The data were better fitted by using the Langmuir model. The Langmuir linear plot of shale is displayed in figure 10b and the obtained parameters for both Langmuir and Freundlich models are reported in Table 5. Isotherms obtained for laterite and sandstone did not display any saturation level which made the prediction of the sorption capacity unreliable. Such a behaviour has already been reported in the literature for phosphate sorption by substrates used in constructed wetlands and phosphate retention in filter materials such as shellsand and filtralite P ® [START_REF] Boeykens | Eutrophication decrease: Phosphate adsorption processes in presence of nitrates[END_REF][START_REF] Drizo | Physicochemical screening of phosphateremoving substrates for use in constructed wetland systems[END_REF]. In any way, the maximum phosphate uptake varied between 0.48 for shale and 8.3 mg g -1 for laterite. Note that the shale dose was ten times higher than those of sandstone and laterite. The ranking order of the maximum sorption capacity was therefore again: laterite  sandstone  shale. Remind that the specific surface area of shale, sandstone and laterite were 6.3, 28.5 and 37.6 m 2 g -1 respectively.

It thus appeared that the material with the lowest specific surface area had the lowest phosphate uptake capacity. One reason for this may be that phosphate sorption was decreased by limited available surface sorption sites in case of shale. Another reason explaining the lowest capacity of shale for phosphate removal was the status of iron. Indeed the shale contained the lowest quantity of iron (Table 1) and the lowest proportion of Fe in iron oxides , i.e. 31 % in comparison with 97 % and 100 % observed for laterite and sandstone, respectively. Iron oxides particles, in particular those with low crystal size, are well known to be preferential adsorption sites for phosphate species. The superiority of laterite for phosphate removal in comparison with shale and sandstone may be correlated to the existence of goethite particles with very low crystal size as testified by the superparamagnetic behaviour observed for this mineral by room temperature Mössbauer spectroscopy. Interstingly, laterite and sandstone were found to exhibit, on a surface area normalized basis, similar sorption capacities close to 0.1 g m -2 (section 2.2.4.), thus suggesting that the difference in reactivity was due to a particle size effect in agreement with Mössbauer analysis.

Phosphate adsorption in hydrodynamic conditions: column tests.

The breaktrough curves of sandstone and laterite obained at a flow rate of 1 mL min -1 and a corresponding HRT of  4 hours are presented in Figure 11. We noticed that the breakthrough curve of laterite presented unsual shape with 2 plateaus. The origin of this behaviour was outside of the scope of the present study. As expected from the batch experiments, the breakthrough observed for laterite occured later than the one observed for sandstone. Indeed, if a threshold of [PO4] = 1 mg L -1 was chosen, i.e. C / C0 of 0.04 with C0 = 25 mg L -1 , the breakthrough occured at V / VP of  147 and 132 for laterite and sandstone, respectively. These values were about 8 times higher than the breakthrough observed previously for shale, i.e. V / VP of  17, even though the mean HRT of ~ 35 hours used for shale was much higher than the one used for sandstone and laterite. Therefore the column experiments confirmed that the capacity of shale to remove phosphate was the lowest as already observed in the batch experiments. The removal capacity at the breakthrough qbreak can be computed to be qbreak = 0.17, 1.23 and 1.30 mg PO4 g -1 for shale, sandstone and laterite respectively. Calcium, magnesium and phosphate ions concentrations are displayed as a function of time (Fig. 12) The initial calcium (4.9 and 9.2 mg L -1 for laterite and sandstone, respectively), magnesium (1.8 and 1 mg L -1 ), aluminum (not detected) and iron (not detected) concentrations in the outflow were in good agreement with the equilibrium data recorded at a pH of 5 in batch experiments. The soluble calcium and magnesium concentrations in the outflow decreased continuously and reached a value lower than 0.5 mg L -1 at the breakthrough. Experiments carried out in hydrodynamic mode required a preliminary step in which the column was pre-equilibrated with ultra-pure water during 24 h before the phosphate solution flowed through the column. The aqueous species concentrations determined at t =0 thus reflected the equilibria solubility that occurred during this period (Fig. 12). In this experimental situation, phosphate removal partially occurred through reaction with the soluble species released in solution during the column equilibration. In fact, a correlation between the decrease in Ca and Mg concentrations and the increase in phosphate concentration at the outflow of the column was evidenced. Iron concentration as a function of time remained below the detection limit (data not shown). Such a behaviour had already been reported in the literature for shellsand and filtralite P  filter materials (Ádám et al., 2007). According to this study, phosphate removal occurred through the accumulation of large amounts of Ca (45-260 mg L -1 ) and Mg ions (5-65 mg L -1 ) at the beginning of the column experiment. The removal capacity at the breakthrough measured for filtralite P  and shellsand can be calculated from the data provided in this previous work (Ádám et al., 2007) and was evaluated to be respectively  0.6 and 0.9 mg PO4 g -1 for a contact time of  4 days. Therefore, the removal capacity of laterite and sandstone from Ivory Coast was at least of the same order of magnitude than the one reported for the shellsand from Hordaland Tørkeri in Norway, which was recognized as a very efficient granular natural material for phosphate removal in hydrodynamic conditions.

Conclusion

The removal of phosphate by three natural materials from Côte d'Ivoire was studied in homogeneous batch suspensions and in static/hydrodynamic conditions in column experiments.

The investigations revealed that all three materials retained phosphate to various extents. Batch experiments showed that the most efficient phosphate adsorbent was laterite followed by sandstone and shale, the sorption efficiency of shale being strongly smaller than those of laterite and sandstone. The removal of phosphate by adsorption on iron oxides was certainly very limited in shale studied here since most of the iron species were situated in structural substitution into the bulk of the solid clay minerals. In contrast most of the iron in laterite and sandstone occurred as in iron oxides. Laterite was also the most efficient mineral for removing phosphate in flow through conditions with a phosphate removal capacity at the breakthrough higher than 1 mg-PO4 g -1 . Therefore, laterite is a very promising mineral for performing extensive dephosphatisation treatments in Ivory Coast and neighbouring countries of Africa [START_REF] Xue | Characteristics and mechanisms of phosphate adsorption onto basic oxygen furnace slag[END_REF].

Figure Captions

Q Q M C C Q Q Q+A A M Q C a.u. Shale M G G Q H Q Q Q a.u. Sandstone H H H Q B G G G Q H Q Q G B Q H G H Q B Q G a.u.
2 theta (degrees)

Fig. 1 -

 1 Fig. 1 -XRD patterns of naturally occurring shale, sandstone and laterite from Ivory Coast. The following abbreviations are used: A: Albite; B: Birnessite; C: Clinochlore; G: goethite; H: Hematite; M: Muscovite and Q: quartz.

Fig. 2 -

 2 Fig.2 -Mössbauer spectra of naturally occurring sandstone (a) and laterite (b) from Ivory Coast at room temperature.

Fig. 3 -

 3 Fig.3 -Mössbauer spectra of naturally occurring shale, sandstone and laterite from Ivory Coast at 12 K (shale) and 9 K (sandstone and laterite).

Fig. 4 -Fig. 5 -Fig. 6 -

 456 Fig. 4 -Phosphate removal efficiency (a) and phosphate uptake (b) for shale, sandstone and laterite as a function of adsorbent dose ([PO4]initial = 25 mg L -1 ; pH ~ 7; time contact: 24 h)

Fig. 7 -

 7 Fig. 7 -Ions concentrations in solution in blank experiments and after reaction with a 25 mg L - 1 initial phosphate solution ([sorbent] = 80 g L -1 ; pH ~ 7; time contact: 24 h).

Fig. 8 -

 8 Fig. 8 -Effect of contact time on phosphate uptake capacity by sandstone, laterite and shale ([PO4]initial = 25 mg L -1 ; pH ~ 7; contact time: 24 h).

Fig. 9 -

 9 Fig.9 -Kinetics modelling of phosphate sorption onto shale (80 g/L) and sandstone and laterite (8 g/L) using pseudo-second order model ([PO4]initial = 25 mg L -1 ; pH ~ 7; time contact: 24 h).

Fig. 11 -

 11 Fig.11-Breakthrough curves obtained for the sandstone and laterite column experiment in hydrodynamic conditions (flow of 1 mL min -1 , C0 = 25 mg PO4 L -1 ).

Figure 12 -

 12 Figure 12 -Ca 2+ , Mg 2+ and PO4 3-concentrations as a function of time in sandstone and laterite

Fig. 1 -Fig. 2 -Fig. 3 -

 123 Fig. 1 -XRD patterns of naturally occurring shale, sandstone and laterite from Ivory Coast. The following abbreviations are used: A: Albite; B: Birnessite; C: Clinochlore; G: goethite; H: Hematite; M: Muscovite and Q: quartz.

Fig. 4 -Fig. 5 -Fig. 6 -Fig. 7 -Fig. 8 -

 45678 Fig. 4 -Phosphate removal efficiency (a) and phosphate uptake (b) for shale, sandstone and laterite as a function of adsorbent dose ([PO4]initial = 25 mg L -1 ; pH ~ 7; time contact: 24 h)

Fig. 9 -Fig 10 -

 910 Fig.9 -Kinetics modelling of phosphate sorption onto shale (80 g/L) and sandstone and laterite (8 g/L) using pseudo-second order model ([PO4]initial = 25 mg L -1 ; pH ~ 7; time contact: 24 h).

Fig. 11 -Fig. 12 -

 1112 Fig. 11 -Breakthrough curves obtained for the sandstone and laterite column experiment in hydrodynamic conditions (flow of 1 mL min -1 , C0 = 25 mg PO4 L -1 ).

Table captions Table 1 -

 captions1 Chemical composition (wt %) of naturally occurring shale, sandstone and laterite samples from Ivory Coast.

Table 2 -

 2 Mössbauer hyperfine parameters of the spectra of shale, sandstone and grey presented in Fig.2. CS: center shift with respect to -Fe at room temperature; : quadrupole splitting or : quadrupole shift ; H hyperfine magnetic field and RA: Relative Area.

Table 3 -

 3 Mössbauer hyperfine parameters of the spectra of shale, sandstone and grey presented in Fig.2. CS: center shift with respect to α -Fe at room temperature ∆: quadrupole splitting or Ԑ: quadrupole shift ; H hyperfine magnetic field and RA: Relative Area.

Table 4 -

 4 Kinetic models parameters for phosphate adsorption on laterite, sandstone and shale after reaction with an initial phosphate concentration of 25 mg L -1 .

Table 5 -

 5 Langmuir and Freundlich isotherm parameters for phosphate adsorption on shale (pH ~ 7; [shale] = 80 g/L; contact time: 24 h)

		SiO2 Al2O3 Fe2O3 MnO MgO CaO Na2O K2O TiO2 P2O5 LoI*
	Shale	56.8 17.5	10.2	0.6	2.1	0.4	0.9	2.4	1	0.1	8.7
	Sandstone 29.3 6.3	54.4	L.D. L.D. 0.1	L.D. 0.1	0.3	1.0	8.7
	Latérite	50.3 14.0	24.5	0.1	0.08	0.08 0.02	0.3	0.7	0.1	9.9
	*Loss on Ignition

Table 1 -

 1 Chemical composition (wt %) of naturally occurring shale, sandstone and laterite samples from Ivory Coast.

	Name	Component	CS	∆ ou Ԑ	H	RA
			(mm.s -1 )	(mm.s -1 )	(kOe)	(%)
	Sandstone	D1 (-FeOOH)	0.32	1.07	-	19
		S1 (-FeOOH)	0.41	-0.17	290	41
		S1 (Fe2O3)	0.37	-0.04	491	40
	Laterite	D1 (-FeOOH)	0.34	0.55	-	55
		S1 (-FeOOH)	0.35	-0.26	340	6
		S1 (Fe2O3)	0.66	-0.1	500	20
		S1 (Fe2O3)	0.36	-0.09	476	19

Table 2 -

 2 Mössbauer hyperfine parameters of the spectra of shale, sandstone and grey presented in Fig.2. CS: center shift with respect to -Fe at room temperature; : quadrupole splitting or : quadrupole shift ; H hyperfine magnetic field and RA: Relative Area.

	Name	T	Component	CS	∆ ou Ԑ	H	RA
		(K)		(mm.s -1 )	(mm.s -1 )	(kOe)	(%)
	Shale	12	D1 (Fe 3+ )	0.470	0.836	-	55
			D2 (Fe 2+ )	1.272	2.71	-	14
			S1 (α-FeOOH)	0.6	-0.2	480	23
			S2 (α -FeOOH)	0.56	-0.08	500	8
	Sandstone	9	S1 (α -FeOOH)	0,48	-0.12	497	51
			S2 (Fe2O3)	0,48	-0.08	532	49
	Laterite	9	S1 (α -FeOOH)	0.46	-0.12	502	28
			S2 (α -FeOOH)	0.46	-0.12	487	32
			S3 (Fe2O3)	0.46	-0.09	532	40
			D1 (Fe(III)))				3

Table 3 -

 3 Mössbauer hyperfine parameters of the spectra of shale, sandstone and grey presented in Fig.2. CS: center shift with respect to α -Fe at room temperature; ∆: quadrupole splitting or

Ԑ: quadrupole shift; H hyperfine magnetic field and RA: Relative Area.

Table 4 -

 4 Kinetic models parameters for phosphate adsorption on laterite, sandstone and shale after reaction with an initial phosphate concentration of 25 mg L -1 .

				Pseudo-first-order kinetic	Pseudo-second-order kinetic
					model			model
	Material	Dose	qexp	qcalc	k1	R 2	qcalc	k2	R 2
		(g L -1 )	(mg g -1 )	(mg g -1 )	(min -1 )		(mg g -1 )	(g mg -1 min -1 )
	Laterite	15	1.70	0.41	0.01	0.703 1.78	0.041	0.999
		8	3.2	1.31	0.01	0.926 3.2	0.006	0.998
	Sandstone	80	0.99	0.29	0.01	0.927 0.99	0.51	1.000
		30	0.87	0.48	0.02	0.963 0.89	0.10	1.000
		8	2.8	1.80	0.01	0.800 2.84	0.01	0.990
	Shale	80	0.21	0.07	0.01	0.955 0.23	0.04	0.999

Table 5 -

 5 Langmuir and Freundlich isotherm parameters for phosphate adsorption on shale (pH ~ 7; [shale] = 80 g/L; contact time: 24 h)

	Model	Equation														Parameters
	Langmuir	e e C q		m e C q		K	q 1 L	m						R 2 qm (mg/g) KL(L/g)	0.995 0.514 0.061
	Freundlich	Ln	q (	e	)		Ln	(	K	f	)		nf 1	Ln	(	C	e	)	R 2 nf	0.955 2.02
																					Kf (L/g)	0,046

Acknowledgements

The authors thank G. Medjahdi (Centre of X-gamma Competence; UMR 7198 CNRS -Université de Lorraine, Nancy-France) for DRX analyses and C. Genois from the LCPME is gratefully acknowledged for ICP analyses.