Rémi Joubaud

Olivier Bernard

Alfred Delville

Alexandre Ern

Benjamin Rotenberg

Pierre Turq

Numerical study of density functional theory with mean spherical approximation for ionic condensation in highly charged confined electrolytes

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Charged solutes and charged interfaces are encountered in numerous natural and technological contexts, largely due to the ability of water to solvate ions and to screen their electrostatic interactions by its high dielectric permittivity. In this paper, we are interested in the study at the nanometric scale of ionic solutions interacting with (or confined by) solid objects carrying a surface charge (clay rocks, colloids, DNA, or electrodes, to name only a few examples). Multi-valent ions and electrolytes near such highly charged objects may depart significantly from the ideal behavior of infinitely diluted solutions.

In bulk situations, the nonideal behavior of ionic solutions arises mainly from two types of effects, which both play a larger role as the ionic concentration increases: long-range electrostatic interactions and short-range excluded volume effects. By treating the former at the mean-field level and ignoring the latter, the pioneering work of Debye and Hückel (DH) [1] identified screening by counterions as a fundamental mechanism for nonideality. The DH theory is valid at low concentrations and high temperature. For more concentrated solutions, further progress has been accomplished within the framework of the primitive models of charged hard spheres in a continuous solvent characterized by its dielectric constant, whereby structural and thermodynamic properties are predicted using integral equation theories, such as the Mean Spherical Approximation (MSA) [START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. I. Method of solution[END_REF][START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. II. Results[END_REF][START_REF] Blum | Mean spherical model for asymmetric electrolytes[END_REF][START_REF] Blum | Mean spherical model for asymmetric electrolytes. 2. Thermodynamic properties and the pair correlation function[END_REF]. Liquid-vapor transition and criticality in bulk ionic solutions have been extensively investigated over the past decades. The properties of the liquid-vapor interface have been studied using a density-functional theory (DFT) within the restricted primitive model (RPM) in which the ions are modelled as equi-sized spheres carrying opposite charges [START_REF] Groh | Liquid-vapor interface of an ionic fluid[END_REF]. Moreover, by combining Bjerrum's ion pairing concept [7] with the MSA or the DH theory through the law of mass action, better descriptions of the critical point and of the liquid-vapor coexistence curve have been obtained [8,9,[START_REF] Fisher | Criticality in ionic fluids: Debye-Hückel theory, Bjerrum, and beyond[END_REF][START_REF] Levin | Criticality in the hard-sphere ionic fluid[END_REF][START_REF] Jiang | Criticality and phase behavior in the restricted-primitive model electrolyte: description of ion association[END_REF].

In the case of confined ionic solutions, electrostatic interactions also control the structure and the phase behavior of the solution [START_REF] Flores-Mena | Criticality in confined ionic fluids[END_REF][START_REF] Pizio | Phase behavior of ionic fluids in slitlike pores: a density functional approach for the restricted primitive model[END_REF][START_REF] Pizio | Phase behavior of the restricted primitive model of ionic fluids with association in slitlike pores. density-functional approach[END_REF]. In some cases, such as clay minerals, cement, ion exchange membranes or lipid bilayers, the counterions compensating the charge of the surface may even be the only ions present in the confined fluid, resulting in a situation similar to the one-component plasma (see [START_REF] Baus | Statistical mechanics of simple Coulomb systems[END_REF] for a review and [START_REF] Jardat | Salt exclusion in charged porous media: a coarse-graining strategy in the case of montmorillonite clays[END_REF] for a discussion closer to the present setting). The starting point for the description of confined ionic solutions is the Poisson-Boltzmann (PB) theory [START_REF] Dresner | Ion exclusion and salt filtering with porous ion-exchange materials[END_REF][START_REF] Dresner | The exclusion of ions from charged microporous solids[END_REF][START_REF] Engström | Ion condensation on planar surfaces. a solution of the Poisson-Boltzmann equation for two parallel charged plates[END_REF]. The PB theory ignores correlations between ions, and no phase transition exists within this approach, although quite high counterion densities can be predicted near charged surfaces in some conditions (e.g., cylindrical geometry [START_REF] Manning | Limiting laws and counterion condensation in polyelectrolyte solutions i. colligative properties[END_REF][START_REF] Bret | Distribution of counterions around a cylindrical polyelectrolyte and manning's condensation theory[END_REF]). For highly charged surfaces or multi-valent ions, a large fraction of the counterions appears condensed near the charged surface, as suggested by Stern to generalize the Gouy-Chapman description of charged surfaces. The remaining ions then feel a much weaker effective charge, which can be described within the PB theory. Nevertheless, the determination of the fraction of condensed ions and the corresponding renormalized charge is not straightforward; results in this direction can be found in [START_REF] Trizac | Alexander's prescription for colloidal charge renormalization[END_REF]. Various approaches have been proposed to incorporate correlations neglected in the PB theory. In the particular case where counterions are the only ions present, a perturbative correction to the PB theory has been established [START_REF] Attard | Beyond Poisson-Boltzmann: images and correlations in the electrical double-layer. I Counterions only[END_REF]. Correlations for confined ions can also be included in integral theories, e.g., within the Anisotropic Hypernetted Chain approximation [START_REF] Kjellander | Attractive double-layer interactions between calcium clay particles[END_REF][START_REF] Kjellander | Double-layer ion correlation forces restrict calcium-clay swelling[END_REF]. Furthermore, the so-called "Strong Coupling" theory allows one to investigate regimes where the interaction with the charged surface is stronger than that between ions [START_REF] Moreira | Binding of similarly charged plates with counterions only[END_REF][START_REF] Netz | Electrostatistics of counter-ions at and between planar charged walls: From Poisson-Boltzmann to the strong-coupling theory[END_REF][START_REF] Yu | The physics of charge inversion in chemical and biological systems[END_REF][START_REF] Samaj | Counterions at highly charged interfaces: From one plate to like-charge attraction[END_REF][START_REF] Samaj | Wigner-crystal formulation of strong-coupling theory for counterions near planar charged interfaces[END_REF] and to explain the origin of the attraction between like-charged surfaces observed under certain conditions [START_REF] Samaj | Counterions at highly charged interfaces: From one plate to like-charge attraction[END_REF][START_REF] Paillusson | Interaction regimes for oppositely charged plates with multivalent counterions[END_REF]. Another successful development for the description of the inhomogeneous primitive model is the use of DFT, which determines structural and thermodynamic properties of an inhomogeneous fluid from a Helmholtz free energy and its functional dependence on the local ionic densities [START_REF] Kralj-Iglic | A simple statistical mechanical approach to the free energy of the electric double layer including the excluded volume effect[END_REF][START_REF] Hansen | Effective interactions between charged colloidal particles : Repulsion, attraction and phase separation[END_REF][START_REF] Netz | Field theory for charged fluids and colloids[END_REF][START_REF] Hansen | Effective interactions between electric double-layers[END_REF][START_REF] Gillespie | Density functional theory of charged, hard-sphere fluids[END_REF][START_REF] Mier-Y Teran | A nonlocal free energy density functional approximation for the electrical double layer[END_REF][START_REF] Löwen | Nonlinear counterion screening in colloidal suspensions[END_REF][START_REF] Oleksy | Wetting and drying scenarios of ionic solutions[END_REF]. Finally, Molecular Dynamics and Monte Carlo simulations have been used to study the properties of bulk and confined electrolytes, described either within the primitive model or with an explicit molecular solvent, thus providing a more realistic description of these complex systems [START_REF] Greberg | Ion-ion correlations in electric double layers from monte carlo simulations and integral equation calculations[END_REF][START_REF] Greberg | Ion-ion correlations in electric double layers from monte carlo simulations and integral equation calculations part 2. case of added salt[END_REF][START_REF] Delville | Electrostatic attraction and/or repulsion between charged colloids: A (NVT) monte-carlo study[END_REF][START_REF] Kjellander | McMillan-Mayer theory for solvent effects in inhomogeneous systems: Calculation of interaction pressure in aqueous electrical double layers[END_REF][START_REF] Pellenq | On the origin of cement cohesion[END_REF]. This paper's goal is to investigate numerically DFT predictions of confined ionic solutions at the nanometric scale by comparing ionic concentration profiles and pressure predictions to those obtained with MC and with the simpler PB approach for the double-layer configuration. We focus here on canonical constraints on the ionic densities fixing their mean value in the fluid domain while ensuring global electro-neutrality. In the present DFT, the Helmholtz free energy of the system depends on the ionic densities and splits into a bulk (ion-ion) and an ion-solid part. The former incorporates the ideal and nonideal effects in the bulk solution regardless of the presence of the charged solid object. The latter is based on a mean-field approach whereby the electrostatic potential is evaluated consistently with the ionic densities as the solution of a Poisson problem with Neumann boundary condition accounting for the surface charges carried by the solid object. Nonideal effects in the bulk free energy account for correlations between ions, including both long-range electrostatic and short-range steric contributions. These correlations are evaluated with the MSA and the Local Density Approximation (LDA) [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]. In particular, the screening length evaluated with the MSA depends on a nondimensional parameter which can be interpreted as a reduced temperature.

⊕ ⊖ ⊖ ⊕ ⊖ ⊕ ⊖ ⊕ ⊖ ⊕ ⊕ ⊖ ⊖ ⊕ ⊕ ⊖ ⊖ ⊕ ⊕ ⊖ ⊕ ⊖ ⊕ ⊕ ⊖ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊖ ⊕ ⊕ ⊖ ⊕ ⊕ L * negatively charged periodic L * Figure 1.
Geometric setting: elementary cell with charged inclusion (left); channel with charged walls (right).

understand the features of DFT-MSA predictions, we perform additional simulations where we consider only the long-ranged electrostatic contribution or the short-ranged steric contribution to the bulk free energy, and we also develop a robust numerical methodology to explore the situation where the reduced temperature falls below the critical temperature predicted by DFT-MSA, leading to phase separation.

The material is organized as follows. In Sect. 2, we present the DFT-MSA considered herein. In Sect. 3, we outline the numerical method for the constrained minimization within DFT and the more direct molecular approach based on MC simulations. In Sect. 4, we compare and analyze predictions obtained with DFT, MC, and PB. In Sect. 5, we present the numerical methodology to treat phase separation within DFT, and we discuss a few selected numerical results. Finally, in Sect. 6, we summarize our main findings and discuss further work.

DFT for confined ionic solutions

We consider ionic solutions in the presence of solid objects carrying surface charges. Two typical settings are that of an ionic solution confined by charged walls and that of an ionic solution in a periodic setting with elementary cells containing a charged inclusion; see Fig. 1. While most of our results focus on the first setting, we briefly consider the more complex geometric setting of charged inclusions to illustrate the capabilities of the present numerical method. The length scale associated with the confinement is denoted by L * and is typically of the order of a nanometer. The domain occupied by the ionic solution is denoted by Ω, and the domain occupied by the solid object is denoted by Ω S . We assume that the solid object carries negative charges with surface density -Σ S (Σ S > 0) on its boundary ∂Ω S .

For simplicity, we consider at most two dissolved ionic species, a counterion (cation) and a co-ion (anion). Ions are modelled within the primitive model as equi-sized spheres of diameter σ. The valences of the ionic species are denoted by Z ± , and the case where Z + + Z -= 0 is referred to as symmetric electrolyte (corresponding to the RPM). The ionic solution is described by the ionic (number) densities ρ = (ρ + , ρ -) in the fluid domain Ω, while the solvent enters the model only by means of its relative dielectric permittivity. In what follows, we consider canonical constraints prescribing the mean-value of the ionic densities in the form

i=± Z i ρ mean i = 1 e|Ω| ∂Ω S Σ S ,
where e is the elementary charge. A simple way of ensuring the global electro-neutrality condition ( 2) is to choose a density of added salt, ρ salt , and to set

(3)

ρ mean + = ρ Σ S + + ρ salt , ρ mean - = -(Z + /Z -)ρ salt ,
where ρ Σ S + = (Z + e|Ω|) -1 ∂Ω S Σ S is the density of neutralizing counterions. In the vanishing limit of added salt, ρ salt → 0, we obtain an ionic solution with only compensating counterions.

The free energy functional of the ionic solution takes the form ( 4)

F(ρ) = F bulk (ρ) + F i-s (ρ).
The bulk free energy functional F bulk accounts for ideal and nonideal ion-ion interactions, while the functional F i-s accounts for ion-solid interactions owing to the surface charges carried by the solid object. The functionals F bulk and F i-s are detailed in Sect. 2.1 and 2.2, respectively.

2.1. Bulk free energy functional. Within the LDA, the bulk free energy is obtained by integrating over the fluid domain Ω the bulk free energy density f , so that (5)

F bulk (ρ) = Ω f (ρ).
The bulk free energy density splits into an ideal part f id and a nonideal part f corr in the form

(6) f (ρ) = f id (ρ) + f corr (ρ), f corr (ρ) = f Coul (ρ) + f HS (ρ),
where f corr accounts for Coulomb electrostatic interactions (through f Coul ) and hard-sphere steric exclusion effects (through f HS ). The ideal part reads

(7) f id (ρ) = 1 βσ 3 i=± ρi (log(ρ i ) -1),
with nondimensional ionic densities ρ± = σ 3 ρ ± , while β = (k B T ) -1 where k B is the Boltzmann constant and T the temperature. The Coulomb contribution is considered to be of the MSA form [START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. I. Method of solution[END_REF][START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. II. Results[END_REF][START_REF] Blum | Mean spherical model for asymmetric electrolytes[END_REF][START_REF] Blum | Mean spherical model for asymmetric electrolytes. 2. Thermodynamic properties and the pair correlation function[END_REF] (8)

f Coul (ρ) = - 1 4πβσ 3 ρα -4Γ MSA (ρ α ) 2 - 16 3 Γ MSA (ρ α ) 3 ,
with ρα = (ρ + /α + ) + (ρ -/α -) and the nondimensional parameters

(9) α ± = σ 4πL B Z 2 ± .
Here, L B stands for the Bjerrum length given by L B = (4πε) -1 βe 2 where ε = ε 0 ε r with ε 0 the vacuum dielectric permittivity and ε r the solvent relative dielectric permittivity. The parameters α ± are related to nondimensional temperatures T * ± such that (10)

T * ± = 4πα ± = k B T σ(4πε) Z 2 ± e 2 ,
while the inverse Debye length κ is such that (11) κσ = ρ1/2 α . The (nondimensional) screening parameter Γ MSA is given by [START_REF] Jiang | Criticality and phase behavior in the restricted-primitive model electrolyte: description of ion association[END_REF] Γ MSA (ρ α ) = 1 2 2ρ 1/2 α + 1

1/2 -1 .
Moreover, the hard-sphere contribution, which hinges on the Carnahan-Starling (CS) expression, is given by

(13) f HS (ρ) =      - 6 πβσ 3 ξ(ρ tot ) 2 (3ξ(ρ tot ) -4) (1 -ξ(ρ tot )) 2 , ξ(ρ tot ) < 1, +∞, ξ(ρ tot ) ≥ 1,
with reduced total ionic density ρtot = ρ+ + ρ-and packing number ξ(ρ tot ) = 1 6 π ρtot . In the dilute limit of low ionic densities, the CS expression can the linearized in the form (CS1) [START_REF] Pizio | Phase behavior of ionic fluids in slitlike pores: a density functional approach for the restricted primitive model[END_REF] f

HS (ρ) = 1 βσ 3 24 π ξ(ρ tot ) 2 .
A simplification occurs for symmetric electrolytes, since, in this case, α + = α -= α and T * + = T * -= T * . Moreover, ρα = α -1 ρtot . As a result, f corr (ρ) only depends on ρtot and its properties can be described in terms of a single reduced temperature T * .

2.2.

Ion-solid free energy functional. The electrostatic interactions between ions and the solid object (as those between ions) are treated at the mean-field level. The electrostatic potential is computed consistently with the ionic densities by solving a Poisson problem in the fluid domain Ω with source term given by the charge density q(ρ) = i=± Z i eρ i and a Neumann boundary condition accounting for the surface charge carried by the solid object on ∂Ω S . Specifically, the electrostatic potential ψ(ρ) solves ( 15)

         -ε∆ψ(ρ) = q(ρ) in Ω, -ε∇ψ(ρ) • n = Σ S on ∂Ω S , ψ(ρ) periodic on ∂Ω \ ∂Ω S , ψ(ρ) Ω = 0,
where n denotes the unit outward normal to Ω. The global electro-neutrality condition [START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. I. Method of solution[END_REF] ensures that problem (15) admits a solution, while the zero mean condition on ψ(ρ) is just one possibility to fix the arbitrariness of ψ(ρ) up to an additive constant. Finally, the functional F i-s (ρ) is given by ( 16)

F i-s (ρ) = 1 2 Ω q(ρ)ψ(ρ) - ∂Ω S Σ S ψ(ρ) = ε 2 Ω |∇ψ(ρ)| 2 ,
the second equality being a direct consequence of ( 15) and integration by parts. Letting G Σ S denote the Green function of [START_REF] Pizio | Phase behavior of the restricted primitive model of ionic fluids with association in slitlike pores. density-functional approach[END_REF] with Dirac function as source term in Ω, the boundary conditions remaining unchanged, the electrostatic potential can be rewritten as a convolution product in the form

(17) ψ(ρ)(x) = Ω G Σ S (x, y)q(ρ)(y) dy.

Numerical methods

This section briefly describes the two numerical methods used in this work, namely the minimization of the free energy described in Sect. 2 and a more direct molecular approach based on Monte Carlo simulations.

3.1. Minimizing the free energy functional. Within the DFT, the equilibrium state of the ionic solution is determined by the following constrained minimization problem: [START_REF] Dresner | Ion exclusion and salt filtering with porous ion-exchange materials[END_REF] ρ minimizes F(ρ) defined by ( 4), ρ satisfies the canonical constraint (1), A (sufficient) well-posedness condition for [START_REF] Dresner | Ion exclusion and salt filtering with porous ion-exchange materials[END_REF] to have a unique solution hinges on the strict convexity of the bulk free energy density f with respect to the ionic densities. This situation, which has been studied mathematically in [START_REF] Ern | Mathematical study of non-ideal electrostatic correlations in equilibrium electrolytes[END_REF] (using the linearized hard-sphere term ( 14)), is encountered when the reduced temperatures T * ± are high enough (typically, above the critical temperature) so that the nonconvexity of the Coulomb interaction term f Coul is compensated by the convexity of the ideal and hard-sphere terms f id and f HS . For a symmetric electrolyte where the function f is univariate, a study of the convexity of this function shows that the critical temperature and (total) density are T * crit ≈ 0.0785 and ρtot,crit ≈ 0.0145 (see Fig. 7 below), in agreement with the result previously derived in [START_REF] Groh | Liquid-vapor interface of an ionic fluid[END_REF]. As is well-known, the MSA overpredicts the critical temperature, which is about T * crit = 0.05 according to Monte Carlo simulations [START_REF] Luijten | University class of criticality in the Restricted Primitive Model electrolyte[END_REF]. For non-symmetric electrolytes, a study of the convexity of the bulk free energy functional can be found in [START_REF] Contento | An efficient algorithm for the double Legendre-Fenchel transform with application to phase separation[END_REF].

The Euler-Lagrange equations associated with [START_REF] Dresner | Ion exclusion and salt filtering with porous ion-exchange materials[END_REF] are formulated in mixed form by introducing as an additional unknown the electrostatic potential ψ = ψ(ρ) solving [START_REF] Pizio | Phase behavior of the restricted primitive model of ionic fluids with association in slitlike pores. density-functional approach[END_REF]. The constraints to be taken into account are the canonical constraint (1) on the ionic densities ρ and the fact that ψ has zero mean value in Ω. In view of finite element discretization, we write the Euler-Lagrange equations in variational form using test functions. Thus, we seek ionic densities ρ, electrostatic potential ψ, and real numbers λ and µ bulk

± such that (19)              Ω ε∇ψ • ∇ϕ + λ ϕ Ω = Ω q(ρ)ϕ - ∂Ω S Σ S ϕ, Ω µ el ± (ρ, ψ)v ± = µ bulk ± |Ω| v Ω , r ψ Ω = 0, s ± ( ρ η,± Ω -ρ bulk ± ) = 0
, for all test functions ϕ and v ± and for all real numbers r and s ± . The electro-chemical potential

µ el ± (ρ, ψ) is such that (20) µ el ± (ρ, ψ) = µ ± (ρ) + Z ± eψ, with chemical potential (21) µ ± (ρ) = ∂ ρ ± f (ρ).
The second equation in [START_REF] Dresner | The exclusion of ions from charged microporous solids[END_REF] expresses the fact that the electro-chemical potentials are constant in Ω (with constant value equal to µ bulk ± ). The Euler-Lagrange equations [START_REF] Dresner | The exclusion of ions from charged microporous solids[END_REF] are discretized using finite elements. The electrostatic potential and the ionic densities are approximated using continuous, piecewise affine functions over a mesh of the fluid domain Ω. This leads to a (large) system of nonlinear algebraic equations which is solved iteratively using a Newton-Raphson algorithm. Each Newton-Raphson iteration consists in assembling a linear system whose solution provides an update for the discrete solution vector. By linearity, at each iteration, the electrostatic potential is evaluated consistently with the ionic densities, and the constraints on the mean value of both the electrostatic potential and the ionic densities are satisfied. A clipping is applied to the update vector resetting to the interval [ǫ, 6/π -ǫ] with small ǫ any ionic density that may have fallen outside this interval during the convergence process, while the ionic densities at the other mesh nodes are rescaled in agreement with the canonical constraint (1). The convergence of the Newton-Raphson algorithm is monitored by checking the Euclidean norm of the discrete residual vector and that of the update vector at each iteration. Typically, convergence is achieved after 15 to 20 iterations, and the is usually needed only at (very) few iterations. More details on the algorithmic aspects can be found in [START_REF] Joubaud | Modélisation mathématique et numérique des fluides à l'échelle nanométrique[END_REF].

3.2. Monte Carlo simulations. Monte Carlo simulations are performed in the Canonical Ensemble to determine the equilibrium distribution of the di-and trivalent ions confined between two charged planes. In each case, the simulation cell is maintained electrically neutral, a condition necessary to reach the thermodynamic limit, i.e., to obtain results independent of the size of the simulation cell. The limiting charged lamellae are structure-less and infinite, with a uniform surface charge density. In the framework of the primitive model, ion-ion and ion-lamella interactions include short-ranged contact repulsion and long-ranged Coulomb potential. The solvent is replaced by a continuum characterized by its relative dielectric permittivity ε r , equal to the relative lamellar dielectric permittivity.

The simulation cell is an asymmetric slit pore, limited in the longitudinal direction by the charged lamellae while periodic boundary conditions are applied in the lateral direction. Since 2D minimum image convention is applied in the transverse directions, the ion-ion, ion-lamella and lamella-lamella electrostatic interactions are cut at the half-width of the simulation cell. In order to reproduce the long-range of the Coulomb potential, an external potential is introduced in a self-consistent manner [START_REF] Van Megen | [END_REF]52,53]. This external electrostatic potential is generated by a supplementary set of infinite planes parallel to the two confining charged lamellae. The uniform surface charge density of these supplementary lamellae is evaluated in a self-consistent manner from the ionic local densities in the simulation cell. A square hole, corresponding to the lateral section of the simulation cell, is cut within each supplementary lamella to remove the contribution from the simulation cell. The convergence of that approximation is carefully checked by doubling the lateral extent of the simulation cell [53]. In the present case, lateral extents of 200 and 400

• A are used successfully for charged interfaces neutralized by divalent and trivalent ions, respectively. Furthermore, three independent and consecutive sets of, at least, 5×10 8 iterations are performed in order to ensure the convergence of the self-consistent procedure. This method has been applied in the literature to investigate the stability of a large class of charged interfaces [54,55,56].

To evaluate the pressure, we proceed as follows. After thermalization of the ionic concentration profiles, the longitudinal component of the pressure tensor is evaluated locally at half separation between the charged limiting surfaces. Such a procedure is valid since, at equilibrium, the longitudinal component of the pressure is constant along the whole interfacial system. In addition to the repulsive entropic contribution evaluated from the local densities of the confined ions, the primitive model includes two additional contributions resulting from the long-range electrostatic couplings and the short-range ion-ion steric repulsions. These two additional contributions result from the ionic correlations leading respectively to antagonistic contributions, i.e., a net attraction for the electrostatic coupling and, obviously, a net repulsion for the ion-ion contact forces. These two components are evaluated by summing the longitudinal component of the net forces transmitted across a fictitious plane surface [55,56].

Results and discussion

In this section we investigate the predictions of DFT-MSA for the ionic density profiles and the pressure, under severe confinement conditions. The results are compared to the reference Monte Carlo (MC) simulations, as well as to the simpler Poisson-Boltzmann (PB) model. In what follows, we consider a confining distance L * = 1 nm (determining the accessible volume) and discuss the influence of the surface charge density Σ S , the ion diameter σ, and salt concentration ρ salt on the ionic profiles and the pressure. We also consider the effect of the relative permittivity ε r and of the valences Z ± for a symmetric electrolyte. In all cases, the reduced temperature T * defined by [START_REF] Fisher | Criticality in ionic fluids: Debye-Hückel theory, Bjerrum, and beyond[END_REF] is in the range [0.08, 0.16], i.e., above the critical point (but not too far from it). In order to understand the origin of the features observed with DFT-MSA, we finally consider separately the effects of the hard sphere and electrostatic correlation contributions to the functional. 4.1. Divalent ions: density profiles. We first examine ionic density profiles for divalent ions (Z ± = ±2), for fixed relative permittivity ε r = 78.5 and salt concentration ρ salt = 0.5 mol L -1 . Fig. 2(a) reports the counter-and co-ions density profiles for a surface charge density Σ S = 0.2 C m -2 and ion diameter σ = 2.55

• A (T * = 0.09). For both profiles, the MC results are better described by PB than by DFT-MSA, in particular near the charged surface where DFT-MSA overestimates the densities. However, the co-ion density is non-monotonous and displays a secondary minimum at the center of the pore. This feature is not reproduced by PB, while DFT-MSA is able to do so. The position of the maximum is nevertheless too close to the surface.

Increasing the surface charge density to Σ S = 0.6 C m -2 , as shown in Fig. 2b, we find that DFT-MSA now underestimates the counterion density at contact. However, DFT-MSA still predicts the observed non-monotonous co-ion profile missed by PB. Comparing the MC co-ion concentration profiles between panels (a) and (b), we observe that the position of the maximum slightly changes, while the peak is more pronounced; the density near the surface and in the center of the pore are reduced accordingly. With DFT-MSA, the position of the peak shifts away from the surface, thereby approaching the MC results. The overestimation of the co-ion density near the surface has now disappeared, and the DFT-MSA prediction is now in very good agreement with the MC results, while a sizable difference remains at the center of the pore.

We now consider, for Σ S = 0.2 C m -2 , an increase in the ion diameter to σ = 4.54

• A (T * = 0.16), as shown in Fig. 2c. While the counterion profile is still better described by PB, DFT-MSA slightly better reproduces the co-ion density. In this regime, the profile increases monotonically from the surface to the center. Compared to PB, DFT-MSA predicts a lower density near the surface and a slightly higher density near the center, i.e., moving in the direction of the MC results. Finally, Fig. 2d reports the results for the same surface charge density and σ = 2.27

• A (T * = 0.08). In this case, the failure of DFT-MSA is dramatic, with a large density of both counter-and co-ions near the surface and accordingly low densities near the center. This behavior is likely due to an incipient phase transition as we approach the critical point. This aspect is further examined in Sect. 5.

4.2.

Divalent ions: pressure. We now turn to the study of the pressure. The pressure from MC calculations is evaluated as described in Sect. 3.2. The pressure from DFT-MSA and PB calculations is evaluated at the center of the pore using the Gibbs-Duhem rules [START_REF] Bret | Distribution of counterions around a cylindrical polyelectrolyte and manning's condensation theory[END_REF] p(ρ

) = ρ + µ + (ρ) + ρ -µ -(ρ) -f (ρ).
For completeness, we also evaluate the pressure in the bulk reservoir in thermodynamic equilibrium with the confined system by evaluating the bulk concentration corresponding to a system at the same chemical potentials µ bulk ± . As the concentration in the bulk have to be equal, we introduce a shift ζ (consistently with our approach for the confined system where we have fixed the mean of the electrostatic potential arbitrarily to 0), so that the two equations satisfied by the densities in the bulk are [START_REF] Trizac | Alexander's prescription for colloidal charge renormalization[END_REF] µ ± (ρ bulk ) + Z ± ζ = µ bulk ± . These equations can be solved by a Newton-Raphson algorithm for each value of the shift parameter ζ, and the actual value of this parameter, which is such that the solution to [START_REF] Trizac | Alexander's prescription for colloidal charge renormalization[END_REF] verifies ρ bulk + = ρ bulk -, is found by dichotomy. Starting from the same system as in Fig. 2a, namely Z ± = ±2, ε r = 78.5, ρ salt = 0.5 mol L -1 , Σ S = 0.2 C m -2 , and σ = 2.55

• A (T * = 0.09), we change T * by increasing σ up to 4.54

• A (T * = 0.16). Results are reported in Fig. 3a. MC simulations indicate that the pressure increases quadratically with T * . PB does not capture the evolution with ion diameter and overestimates the pressure in this range. In contrast, DFT-MSA accurately describes the change in pressure with T * even though its predictions are shifted (by an almost constant value in this range) from the MC results. The PB and DFT-MSA pressures for bulk electrolytes in thermodynamic equilibrium with the confined system are also displayed and found to be almost identical to the pressures under confinement. The effect of the surface charge density can be seen in Fig. 3b. Increasing Σ S to 0.4 C m -2 , we observe that the MC pressure decreases at low T * , down to negative values for T * ≤ 0.11, and increases for the largest values of T * . While the constant PB result increases with Σ S , DFT-MSA better reproduces the quadratic variation predicted by MC simulations.

In Fig. 4, T * is changed by increasing the relative permittivity ε r for a fixed ion diameter σ. The pressure increases almost linearly with T * . The PB pressure now also increases, but with too small a slope and starting from a much larger value. Instead, DFT-MSA is almost in quantitative agreement with MC for a salt concentration ρ salt = 0.5 mol L -1 . Here again, the pressure in the confined system is almost identical to that in the bulk reservoir with which it is in equilibrium. In the absence of salt, the pressure decreases for MC, DFT-MSA, and PB. MC simulations indicate that it even becomes negative for T * ≤ 0.11, a feature that is not reproduced by either PB or DFT-MSA. The deviations of DFT-MSA from MC results are much larger than in the presence of salt, but decrease as T * increases.

Trivalent ions.

We consider now trivalent ions (Z ± = ±3) at the same ionic strength as for the above divalent ions, leading to a salt concentration ρ salt = 0.2222 mol L -1 , in a solvent with relative permittivity ε r = 78.5 and ion diameter σ = 5.1

• A (T * ≈ 0.08 (slightly) above the critical point). The surface charge density is Σ S = 0.2 C m -2 . The ionic density profiles are reported in Fig. 5a. The MC results are qualitatively similar to the ones for divalent ions in this regime (see Fig. 2c). The predictions of DFT-MSA are worse than those of PB for both counter-and co-ions, except for the stronger exclusion of co-ions very close to the surface. The pressure, reported in Fig. 5b as a function of T * (changed by increasing ε r ), behaves as in the divalent case. Despite its poorer prediction of the ionic profiles, DFT-MSA reproduces much better than PB the pressure predicted by MC. 4.4. Hard sphere and electrostatic contributions. We finally investigate the origin of the features observed on the ionic density profiles within DFT-MSA. Compared to PB, the MSA free energy functional includes two additional contributions: the hard sphere (hs) volume exclusion term defined in (8) and the electrostatic correlation (el) term defined in [START_REF] Flores-Mena | Criticality in confined ionic fluids[END_REF]. In order to assess their effects on the full MSA solution, we consider two additional cases, where we include only one or the other term in the functional. Results for counter-and co-ions are reported in Figs. 6a and6b, respectively, for divalent ions (Z ± = ±2), Σ S = 0.2 C m -2 , ρ salt = 0.5 mol L -1 , σ = 2.55

• A, and ε r = 78.5 (T * = 0.09), i.e., the same situation as Fig. 2a. Adding only the hard sphere term to PB has almost no effect on the co-ions, but slightly decreases the counterion concentration near the surface, as expected, thereby providing a better agreement with the MC results. On the contrary, adding only the electrostatic correlations (at the MSA level) results under the present conditions in a wild overestimate of both densities near the surface. In the absence of hard core repulsion, the chosen T * probably corresponds in the bulk to a phase coexistence, and the separation between a high density (close to the surface) and a low density (near the center) region is reminiscent of this transition. Even though this pure electrostatic correction is not relevant to the description of the real system, it suggests that the full DFT-MSA overestimate of the counterion density near the surface can be traced back to this contribution. In addition, these electrostatic correlations are at the origin of the non-monotonous behavior of the co-ion density (when it occurs). Finally, since the hard sphere contribution can only increase the pressure, the electrostatic correlation term is responsible for the smaller pressure with DFT-MSA than with PB, resulting in the better agreement with MC results.

Phase separation within DFT-MSA

This section briefly examines the present DFT-MSA when the bulk free energy density is no longer convex with respect to the ionic densities. We recall that in the case of bulk ionic solutions, this regime leads to liquid-vapor transition. We emphasize that the present DFT-MSA results substantially differ from MC predictions. 5.1. The constrained minimization problem. In confined ionic solutions, the presence of the charged solid object induces gradients in the ionic densities and, thereby, departures from local electro-neutrality condition i=± Z i eρ i = 0. As a result, the whole state space K := {ρ ∈ R 2 ; ρ ± ≥ 0; ξ(ρ) < 1} can potentially be explored by the ionic densities at different points in Ω. Letting f * * be the convex hull of f , we introduce the subset ( 24)

K * * = {ρ ∈ K; f (ρ) = f * * (ρ)},
and we supplement the constrained minimization problem [START_REF] Dresner | Ion exclusion and salt filtering with porous ion-exchange materials[END_REF] with the additional constraint that, everywhere in Ω, the ionic densities take values in the subset K * * . This constraint hinges on the assumption that the bulk system should be thermodynamically stable regardless of confinement; this is why the convexification is performed only on the bulk part of the free energy functional. The constrained minimization problem can be stated as follows: The mathematical analysis of [START_REF] Kjellander | Attractive double-layer interactions between calcium clay particles[END_REF] goes beyond the present scope. There is an extensive mathematical (and numerical) bibliography on phase separation in other settings than confined electrolytes. In general, the functional used in such phase-field theory is a double-well potential inducing a partition of the state space into two phases. Most of the studies have been performed in bulk situations where the state variable (here, the ionic densities) is constant in each phase. The key idea is to regularize the problem by minimizing a regularized functional where a least-squares penalty on the gradient of the state variable is added to the double-well potential. For the mathematical analysis showing that the limit minimization problem as the regularization parameter tends to zero involves the convex hull of the double-well potential and that this limit represents a choice criterion among minimizers of the original nonconvex variational problem so as to minimize interfaces by sampling states where the potential and its convex hull coincide, we refer, e.g., to the seminal work in [START_REF] Modica | The gradient theory of phase transitions and the minimal interface criterion[END_REF]; further references and links to Γ-convergence can be found in [START_REF] Braides | Γ-convergence for Beginners[END_REF]. The situation where the double-well potential is perturbed by, e.g., a nonlocal operator is more delicate to analyze; we mention [START_REF] Muratov | Droplet phases in non-local Ginzburg-Landau models with Coulomb repulsion in two dimensions[END_REF][START_REF] Choksi | Periodic phase separation: the periodic isoperimetric and Cahn-Hilliard problems[END_REF], and, more recently, [START_REF] Goldman | The Gamma-limit of the two-dimensional Ohta-Kawasaki energy. I. Droplet density[END_REF] where a nonlocal operator accounting for Coulomb interaction is added.

In the case of symmetric electrolytes, it is possible to approximate with very good accuracy the subset K * * . To this purpose, we write f (ρ) = f id (ρ) + f corr (ρ tot ) exploiting the fact that f corr only depends on ρ tot . Re-arranging the ideal part leads to f (ρ) = f tot (ρ tot ) + f rel (ρ), where [START_REF] Kjellander | Double-layer ion correlation forces restrict calcium-clay swelling[END_REF] f tot (ρ tot ) = 1 βσ 3 ρtot (log(ρ tot ) -1) + f corr (ρ tot ), and [START_REF] Moreira | Binding of similarly charged plates with counterions only[END_REF] f rel (ρ) = 1 βσ 3 i=± ρi log ρi ρtot .

The bivariate function f rel is convex in the ionic densities, and the subset K * * is approximated by studying the convexity properties of the univariate function f tot (i.e., we assume that f * * ∼ f rel + f * * tot ). The convexity properties of f tot are deduced from the phase diagram of Fig. 7: for a given reduced temperature T * below the critical temperature T * crit , the binodal points of f tot , denoted by ρ♭ tot and ρ♯ tot , are determined from the liquid-vapor coexistence curves of Fig. 7. The approximate subset K * * of the two-dimensional phase space then consists of two connected components (associated with the dilute and condensed phases) separated by the stripe {ρ ♭ tot ≤ ρ+ + ρ-≤ ρ♯ tot } whose boundary are two parallel lines in the (ρ + , ρ -)plane. A comparison with the exact subset K * * determined numerically by computing the double Legendre-Fenchel transform of f (see [START_REF] Contento | An efficient algorithm for the double Legendre-Fenchel transform with application to phase separation[END_REF]) shows relative differences in predicted ionic densities below 10 -4 away from the interface and 10 -2 near the interface. 5.2. Numerical method. Assuming (25) to be well-posed, we devise a numerical method to approximate its solution. Following the above discussion, we introduce the following regularization of the free energy functional:

(28) F η (ρ) = F(ρ) + η 2 L 5 * β i=± Ω 1 2 |∇ρ i | 2 ,
with nondimensional regularization parameter η > 0. The constrained minimization problem ( 25) is then replaced by seeking ionic densities ρ η such that (29) ρ η minimizes F η (ρ η ) defined by [START_REF] Netz | Electrostatistics of counter-ions at and between planar charged walls: From Poisson-Boltzmann to the strong-coupling theory[END_REF],

ρ η satisfies the canonical constraint (1).
This problem is solved numerically using the techniques outlined in Sect. 3.1, where, in the Newton-Raphson iterations, the chemical potentials are evaluated using the convex hull f * * of the bulk free energy density. We assume that, similarly to the simpler setting analyzed in [START_REF] Braides | Γ-convergence for Beginners[END_REF], the singularly perturbed limit η → 0 + represents a choice criterion between all minimizers of the nonconvex minimization problem (18) by sampling thermodynamically stable states where f and f * * coincide. Solving [START_REF] Yu | The physics of charge inversion in chemical and biological systems[END_REF] with a fixed positive value of the regularization parameter η typically leads to a diffuse interface between the two phases, that is, the set I := {x ∈ Ω; ρ(x) ∈ K * * } has positive measure. The sharp interface approximation η → 0 + leads to two distinct phases occupying the whole fluid domain Ω, while the measure of I tends to zero. The parameter η plays herein the role of a numerical parameter; from a physical viewpoint, a diffuse interface with a fixed value of η can be more appropriate. We refer to [START_REF] Joubaud | Modélisation mathématique et numérique des fluides à l'échelle nanométrique[END_REF] for further insight concerning the numerical methodology. 5.3. Results and discussion. Fig. 8 presents concentration profiles of counterions and coions as a function of the distance to the wall in the double-layer configuration (L * = 1 nm). We consider a symmetric divalent (left panel) or trivalent (right panel) electrolyte with Σ S = 0.2 C m -2 , 0.4 C m -2 , and 0.6 C m -2 , ρ salt = 0.5 mol L -1 , and σ = 2.27

• A (T * = 0.0775). We observe that a single phase transition occurs within the pore, and that the position of the interface shifts toward the center as the surface charge density is increased. In all cases, the condensed phase contains most of the counterions and co-ions. In the divalent case with moderate ion diameter, the co-ion profile is non-monotone in the condensed phase, a feature that almost disappears in the trivalent case with large diameter. Fig. 9 presents the cloud of points (ρ + , ρ -) in the state space for the same parameters as in Fig. 8a. Each point corresponds to the value obtained at a node of the finite element mesh. For fixed Σ S , moving from right to left along a set of points corresponds to moving away from the charged wall. The two parallel oblique lines indicate the boundary of the subset K * * in which the ionic densities are sought. We notice that the constraint that the ionic densities take values in the subset K * * is very well satisfied by the present numerical method, up to very few mesh nodes (less than 0.1%). Finally, to illustrate the capabilities of the numerical method to handle more complex geometries, Fig. 10 presents isocontours of counterion and co-ion concentrations for the configuration with circular inclusions with L * = 1 nm, Σ S = 0.2 C m -2 , and ρ salt = 0.5 mol L -1 . We observe again the formation of a single interface whose shape closely follows that of the circular inclusion, and that the condensed phase contains most of the counterions and co-ions.

Conclusions

The main findings of the present study can be summarized as follows. Ionic concentration profiles are generally better described by PB than by DFT-MSA under the conditions investigated here, namely multi-valent ions under strong confinement, although DFT-MSA is able to capture the non-monotonicity of co-ion profiles missed by PB. The relatively accurate predictions obtained with PB in some situations can be partly attributed to the presence of added salt which tempers the effects of the confining double-layer. Furthermore, the failure of DFT-MSA for the lowest reduced temperatures considered herein is related to an incipient phase separation predicted by DFT-MSA at too high a critical temperature in the bulk. The conclusions regarding the pressure are quite different. Pressure predictions with DFT-MSA are more accurate than with PB, since DFT-MSA reproduces pressure variations with the reduced temperature and, at the same time, yields values closer to those predicted by MC. These results confirm that nonideal effects are important to describe highly confined ionic solutions in thermodynamic equilibrium with a reservoir.

More elaborate models for nonideality can be incorporated into the present DFT in future work. Improved free energy functionals can be handled within the present numerical methodology if they are known in closed analytic form or in tabulated form on a sufficiently fine set of values for the ionic densities allowing for reasonably accurate data interpolation. These can include improvements in the description of electrostatic correlations near the critical point, of excluded volume effects near the solid surface (which can lead to attraction between likecharged surfaces), and of ion-solvent interactions especially for multi-valent ions and very small confinements. In isotropic media, accounting for Bjerrum's ion pairing improves the description of the coexistence curve and the value of the critical point. To better describe the strong electrostatic interactions between the ions and charged surfaces at low temperature, a formulation of the strong coupling theory has been obtained from the Wigner crystal [START_REF] Samaj | Wigner-crystal formulation of strong-coupling theory for counterions near planar charged interfaces[END_REF]. An interpolation at short distance between the strong coupling and Poisson-Boltzmann theories should provide a better account of the correlations in the condensed ion profiles. Another interesting approach is to resort to a Weighted Density Approximation (WDA) instead of the LDA whereby a coarse-grained density obtained by convolution of the local density is used in the free energy; see, e.g., [START_REF] Oleksy | Wetting and drying scenarios of ionic solutions[END_REF].

The present methodology, which is applicable to multi-dimensional geometries, can be of interest to study counterion condensation in confined ionic solutions in other applicative settings. As natural extensions of this work, one could now include the effect of nonideality in multi-scale approaches which currently use the PB theory as a starting point [START_REF] Moyne | A two-scale model for coupled electro-chemo-mechanical phenomena and onsager's reciprocity relations in expansive clays: I homogenization analysis[END_REF][START_REF] Moyne | A two-scale model for coupled electro-chemo-mechanical phenomena and onsager's reciprocity relations in expansive clays: II computational validation[END_REF][START_REF] Allaire | Homogenization of the linearized ionic transport equations in rigid periodic porous media[END_REF][START_REF] Allaire | Role of non-ideality for the ion transport in porous media: derivation of the macroscopic equations using upscaling[END_REF]. Ionic fluids in the presence of random confinement can also be studied by sampling realizations of the geometric configuration of the nanopores and the results compared to previous work, e.g., [START_REF] Pellicane | Replica Ornstein-Zernike self-consistent theory for mixtures in random pores[END_REF]. Another important issue is that of dynamical properties, in particular the coupling of ion dynamics with fluid flows resulting in electro-kinetic effects [START_REF] Pagonabarraga | Recent advances in the modelling and simulation of electrokinetic effects: bridging the gap between atomistic and macroscopic descriptions[END_REF]. The present DFT formulation can be used to include nonideality in numerical approaches which couple hydrodynamic flows to ionic fluxes via time-dependent DFT [START_REF] Rotenberg | Coarse-grained simulations of charge, current and flow in heterogeneous media[END_REF].

( 1 )

 1 ρ ± Ω = ρ mean ± , where • Ω denotes the mean-value of a function in the fluid domain Ω. The prescribed quantities ρ mean ± must satisfy the global electro-neutrality condition (2)

Figure 2 . 54 •

 254 Figure 2. (Color online) Density profiles for divalent ions (Z ± = ±2), relative permittivity ε r = 78.5, and salt concentration ρ salt = 0.5 mol L -1 . (a) Surface charge density Σ S = 0.2 C m -2 and ion diameter σ = 2.55 • A (T * = 0.09). (b) Σ S = 0.6 C m -2 and σ = 2.55 • A (T * = 0.09). (c) Σ S = 0.2 C m -2 and σ = 4.54 • A (T * = 0.16). (d) Σ S = 0.2 C m -2 and σ = 2.27 • A (T * = 0.08). Solid lines and + indicate counterion profiles, while dashed lines and × correspond to co-ions. DFT-MSA (blue) is compared to Poisson-Boltzmann (black) and Monte Carlo (symbols) results.

Figure 3 .

 3 Figure 3. (Color online) Pressure for divalent ions (Z ± = ±2), relative permittivity ε r = 78.5, and salt concentration ρ salt = 0.5 mol L -1 , as a function of reduced temperature T * changed by increasing the ion diameter σ from left to right, for a surface charge density Σ S = 0.2 (a) and 0.4 C m -2 (b). DFT-MSA (blue) is compared to Poisson-Boltzmann (black) and Monte Carlo (symbols) results. Dashed lines indicate the DFT-MSA and PB results in the bulk in thermodynamic equilibrium with the system.

Figure 4 .

 4 Figure 4. (Color online) Pressure for divalent ions (Z ± = ±2), surface charge density Σ S = 0.2 C m -2 , and ion diameter σ = 2.27 • A, as a function of reduced temperature T * changed by increasing the relative permittivity ε r from left to right, for salt concentration ρ salt = 0.5 (a) and 0 mol L -1 (b). DFT-MSA (blue) is compared to Poisson-Boltzmann (black) and Monte Carlo (symbols) results. Dashed lines indicate the DFT-MSA and PB results in the bulk in thermodynamic equilibrium with the system.

ρFigure 5 . 1 •

 51 Figure 5. (Color online) Trivalent ions (Z ± = ±3) with ion diameter σ = 5.1 • A, salt concentration ρ salt = 0.2222 mol L -1 and surface charge density Σ S = 0.2 C m -2 . (a) Ionic density profiles for a relative permittivity ε r = 78.5 (T * ≈ 0.08); see Fig. 2 for legends. (b) Pressure as a function of reduced temperature T * changed by increasing the relative permittivity ε r from left to right; see Fig. 3 for legends.

Figure 6 . 55 •

 655 Figure 6. (Color online) Counterion (a) and co-ion (b) density profiles for divalent ions (Z ± = ±2), surface charge density Σ S = 0.2 C m -2 , salt concentration ρ salt = 0.5 mol L -1 , ion diameter σ = 2.55 • A, and relative permittivity ε r = 78.5 (T * = 0.09). Monte Carlo results (red ×) are compared to PB (dotted black line), full DFT-MSA (solid blue line) and DFT with only the hard-sphere contribution (hs, black dashed line).

Figure 7 .

 7 Figure 7. Phase diagram for a bulk symmetric electrolyte. For each T * below the critical temperature, the curves provide the two binodal values for the reduced total density.

Figure 8 .

 8 Figure 8. (Color online) Concentration of counterions (black) and co-ions (blue) as a function of the distance to the wall for the double-layer configuration (L * = 1 nm) with Σ S = 0.2 C m -2 , 0.4 C m -2 , and 0.6 C m -2 ; (a) divalent electrolyte with ρ salt = 0.5 mol L -1 and σ = 2.27 • A (T * = 0.0775) and (b) trivalent electrolyte with ρ salt = 0.2222 mol L -1 and σ = 4.5 • A (T * = 0.0705).

Figure 9 .

 9 Figure 9. Double-layer configuration (L * = 1 nm) for a symmetric divalent electrolyte with Σ S = 0.2 and 0.6 C m -2 , and ρ salt = 0.5 mol L -1 . Each point corresponds to ionic densities obtained at a certain node of the finite element mesh. The two parallel oblique lines indicate the boundary of the subset K * * .

Figure 10 .

 10 Figure 10. Circular inclusion within symmetric divalent electrolyte for L * = 1 nm, Σ S = 0.2 C m -2 , and ρ salt = 0.5 mol L -1 ; concentration isocontours for (a) counterions (values between 0.27 and 25.3 mol L -1 ) and (b) co-ions (values between 0.09 and 1.82 mol L -1 ).

Acknowledgments

Partial support by ANDRA (French National Agency for the Management of Nuclear Waste) is gratefully acknowledged; in particular, the first author has been supported through a PhD fellowship by ANDRA.


Q Zhang

I Deniaud

E Caillaud

C Baron

DESCRIPTIVE MODEL FOR INTERPRETING INNOVATIVE DESIGN

Keywords: innovative design, descriptive model, system 1Introduction

. So it is difficult to clarify design problems at the beginning of design. Moreover, these traditional models cannot make a clear distinction between routine design and innovative design. That is, there is no descriptive model to clearly reflect the creative activities. It was alsorecognized that these traditional representations do not adequately represent either the designer or the environment. These failings result in part from the assumption that the designer is just apassive executor, instead of an active one. Design is a process during which the designer performs actions in order to change the environment. So the designer observes and interprets what they are "seeing", and then decide on new actions, instead of simply performing the recommended actions. In order to address these limitations, we introduce a descriptive model of innovative design. The model will reflect the actual process and pattern of innovative design,will locate innovation opportunities in the design process and support a systematic perspective whose focus is on the external and internal factors affecting the success of innovative design. In section 2, we discuss the background of the descriptive model by defining innovative design and treating design process as a system. In section 3 we present our proposal for a descriptive model which is organized to achieve the stated goals. In section 4 we presentour concluding remarks.

Background

What is innovative design?

Defining innovative design necessarily requires an understanding of the concepts of innovation and design. Design has not been adequately defined: researchers agree that it is a process, but disagree on what kind of process it is. Some have considered it as a rational problem-solving process, others as a reflective process, and still others as an evolving process between knowledge and concept. There is more agreement about innovation, namely, that it consists of two components: creativity and implementation. Creativity refers to how to generate more and better ideas; implementation concerns how to translate the ideas into product. Looking at one of the definitions of innovative design, one could define it with design variables and forms. However, this still does not go beyond the known framework. Another definition expresses innovative design as a means to break away from the existing rules and generate new rules. The former just class innovative design into one of design categories from variable and form, and does not view innovation as the core part of it. The latter enlarges the scope of innovative design, and refers to a new form of design combining the traditional design function with new elements. However, a new form of design involves not only the product, but also organization and process. We can define innovative design, then, as some kind of process that applies a creative idea to create a product, process, or service for a customer and market. An innovative design should break away existing forms, and demonstrate three characteristics: (1) Novelty. The result of innovative design is different from all previously existing products.

(2) Value. The value of innovative design is related to a human purpose, and should be judged by the customer and society. (3) Commercialization. Innovative design is distinguished from the term creative design because it involves commercial transactions.

Treating the process of innovative design as a system

Before moving further, we need a better understanding of the term "process" and "system". A process is "a set of interrelated or interacting activities which transform inputs into outputs" (IEEE Std 15288 2004). Any work done to produce a result, including creativity and innovation, involves a process. Thus, we can think of the work on innovative design as a process. Without understanding the underlying process, it is difficult to know how the results were achieved, or whether they were good or bad. A system is "a combination of interacting elements organized to achieve one more stated purposes" (INCOSE 2006), or "a set or arrangement of elements that are related, and whose behaviour satisfies operational needs and provides for the life cycle sustainment of the products" (IEEE Std 12202005). Comparing the definition of a process with the definition of a system, we can conclude that a process is a kind of system. Thus, the process of innovative design can be viewed as a system. The process model of innovative design is a typical system model. The literature on the system model is extensive, although only a small subset of this literature focuses on innovative design. However, it is necessary to summarize some fundamental propositions which describe the modelling process of innovative design based on these models.

 Design is something of an art but with many consistent patterns Similarly, while innovative design seeks to design something innovative, the designer or the design team tends to follow a pattern. That is, the process of innovative design has some repeatable structures.  The process model of innovative design is to plan and control the process, and determines design activities. Structured methods are valuable for managing these design activities.  The innovation cannot be fully mechanized. Despite the advances in understanding of innovation, artificial intelligence and techniques, we still cannot accurately prescribe a comprehensive mechanism for innovation. Although these propositions cannot ensure our success in the modelling process, they provide us with the theoretical basis for a descriptive model of the innovative process.

Descriptive model

Based on the propositions in Section 2, in designing the model, we have had three goals in mind.

(1) The model should represent the fundamental processes of innovative design to be useful in understanding the patterns involved in the process.

(2) The model should involve all the necessary elements in the innovative design process, so that we can know when and where the innovations could emerge, proceed and grow.

(3) The model should be comprehensive enough that it will be useful to the designers when discussing different factors affecting the success of innovative design. We believe that the model addresses these goals reasonably well and can, therefore, be useful in interpreting innovative design. Therefore, according to the three goals, the main question of a descriptive model becomes "How should we describe the process of innovative design." This question can be divided into two sub-questions: "What processes should be performed", and "How should we structure these processes." As for the first sub-question, we construct the basic framework of innovative design based on the FBS(Function-Behaviour-Structure) model (Section 3.1). However, the basic framework is not able to take into consideration the entire innovative design process. Referring to the chained-linked model of innovation,the basic framework is extended by considering the innovation process(Section 3.2.1). Because the structure of innovative design is not linear, and inspired by the model of the CPS(creative problem process), we solve the second sub-question by placing the basic framework into a circular form (Section 3.2.2). Further, we introduce the designer and the environment space into the model (Section 3.3). Section 3.4 integrates the two aspects, and finally proposes the descriptive model of innovative design.

Constructing the basic framework of innovative design

John Gero proposed his FBS model of design as a theoretical base for understanding design, and as a conceptual basis for computerized tools intended for supporting the designer [START_REF] Gero | Design prototypes: a knowledge representation schema for design[END_REF]]. According to the FBS model, design consists of eight elementary steps, which are defined in terms of the key concepts of function, behaviour and structure. Five of them transform the desired functions of the artefact into the design descriptions in sequence. And the designer carries out these steps on the basis of knowledge stored in and retrieved by design prototypes. Vermaas and Dorst [START_REF] Vermaas | On the conceptual framework of John Gero's FBS-model and the prescriptive aims of design methodology[END_REF]] identify two problems with the FBS model: the absence of a stable definition of function, and the model"s double aims of describing the actual design and prescribing improved design. Because the second problem involves explaining how descriptive and prescriptive modelling could be distinguished and connected,it goes beyond the scope of this paper. So we focus on a solution to treat the first problem. Vermaas and Dorst"s solution involves a redefinition of the three basic concepts (Table 1), and a reformulation of the eight elementary processes (Table 2). These proposals and analyses are considered as useful ones for understanding and developing the FBS model [START_REF] Per | The ontology of Gero's FBS model of designing[END_REF]], particularly in (1) ensuring the stable conceptual framework over time, and (2) a clear separation between the internal concept (purpose) and the non-internal concept (function, behaviour and structure). The new FBS model modified by Vermaas and Dorst has been chosen as constructing the basic framework of innovative design, because:

1. The new FBS model can be very useful in integrating the design process and the creative process [START_REF] Howard | Describing the creative design process by the integration of engineering design and cognitive psychology literature[END_REF]]. The creative process is an integral part of the innovative design. Three creative process elements (analysis, generation and evaluation) can be mapped onto a view of the design process (Table 2). Tables 2 shows that Processes 1, 2, 6, 7, 8 correspond with the generation stage of the creative process, which results in the generation of new function, new behaviour or new structure. The representation linking the creative process and design process provides us with a way, by which it is possible to judge whether the design stage properly corresponds to an innovative process. It will enable us to focus on the creative activities that the designer is actually performing. 2. The new FBS model accords with the system view of innovative design in Section 2.2, and it appears to be very appropriate to describe the main characteristics of innovative design. A foundational concept of the FBS model is the concept of "design state". A design state is defined as a representation of all possible states which also can be seen as a representation of all possible solutions at a certain moment in time. Because the design states are transferred with time, the state-transition reflects directly the evolving process between design problem and design solution in innovative design. Moreover, the transitions of these design states include a series of design activities. 3. The degree of innovation depends on the distance from the existing design output with the variables and the goals of design, the new FBS model clearly represents these parameters.

According to the definitions in Table 1, the "variables" of design can be assigned to the "behaviour", the "structure" and the "function", while the "goal" of design corresponds to the "purpose". The greater the difference between these parameters, the higher the innovation degree.

Table 1 The definitions of the key concepts on Vermaas and Dorst's view

Purpose(P)

The design intention Function(F) Those physical dispositions of an artefact that contribute to the purpose for which the artefact is designed Behaviour(B)

The physical disposition of the artefact

Structure(S)

The materials of the artefact, the dimensions and geometry of these materials, and their topological relations

Table 2 The eight elementary processes of the FBS model

The eight elementary process [START_REF] Vermaas | On the conceptual framework of John Gero's FBS-model and the prescriptive aims of design methodology[END_REF] Activities in the creative process [START_REF] Howard | Describing the creative design process by the integration of engineering design and cognitive psychology literature[END_REF] Therefore, the new FBS model can provide us with a framework to describe the process of innovative design. Based on the definitions and the processes of Tables 1 and2, the basic framework for innovative design is constructed, as showed in Figure 1. Figure1. The original basic framework of innovative design based on the FBS model Although the new FBS model provides us with a framework, it lacks relevant elements necessary to completely describe the entire process system. Reflection is an important concept for design and innovation, by which the designer can interact with the internal and external environments and when appropriate propose new or even innovative solutions. The sources of information for innovation include internal factors, external/commercial factors, education /research establishments, generally available information and other external factors. However, in its current state, the new FBS model can simply focus on the internal innovation of design process, but cannot consider other possible innovation opportunities associated with the external. The next two sections introduce the relevant elements (namely innovation, the designer and the environment space) which are not taken into account in the FBS model, but these elements are essential and necessary. After defining these elements, we consider how to incorporate them into the descriptive model we propose.

Modifying the basic framework by considering the innovation process

As we have discussed in Section 2.1, an innovative design should break away from existing forms, and demonstrate three characteristics: Novelty, Value and Commercialization. These characteristics are relevant with innovation process. Therefore, the importance of an understanding of how to combine innovation process with design shapes the ways in which we model innovative design. The innovation process concerns why and how innovation emerges, proceeds and grows. There is not a widely agreed upon model of a firm-level process of innovation that incorporates all the requirement of different disciplines. In practice, innovation is a coupling and matching process, and involves many iterations and much feedback. Numerous recent works recognize these characteristics, and propose more complex, variegated and interactive models. In the section, we will modify the basic framework to adapt to these characteristics.

Extending the basic framework by introducing the requirement

One alternative to the linear model, called the "chained-linked model", is proposed by [START_REF] Kline | The Positive sum strategy: harnessing technology for economic growth[END_REF]. They describe how feedback and loops allow the potential innovator to seek existing inter-firm and intra-firm knowledge in the innovation process. Figure2 is a redrawn version of the chained-linked model. The chained-linked model presents a great advance in explaining the necessary interaction for successful innovation. According to the model, a key element in determining the success of innovation is the extent to which a firm manages to maintain effective links between these phases of innovative process. As for innovative design, these effective links and feedback loops provide the probablesources of innovation. These could include, for example, the feedback from the market, technical advances, and the accumulated knowledge. Recall that the basic framework of innovative design that was shown in Figure 1. The original framework, due to the lack of a comprehensive representation of the sources of innovation, is restricted to the internal process of innovative design. According to the descriptions of these elementary processes, these creative processes (Processes 1, 2, 6, 7, 8) result in the innovative outputs, namely new functions (by Processes 1, 7), new purposes (by Process 8) and new structures (by Processes 2, 6). In terms of an innovation process, these creative processes that produce the innovative outputs are not comprehensive. In light of the chained-linked model, the new purposes do not only result from the process S→P(Process 8), but also from technical push and market pull. Technical push describes the scientific and technological advances that push a new idea into the design process, while market pull places more emphasis on market needs.

To better and comprehensively exploit this process, we add requirements (R) into the original framework of innovative design. First, the designer analyzes the requirements(R), which come from the technical advance and the market needs, leading to the interpreted purposes (P). Then the purposes (P) are augmented by the designer"s own experience or the firm"s accumulated knowledge. Therefore, the new basic framework of innovative design is shown in Figure 3. An initial design state space is built up by the requirements (R). Because the transformation from R to P (Process 0) brings on a new purpose, the process 0 is also the generation stage of the creative process.

Circularity in the basic framework

In fact, the new version is still a logical development of the innovative design shown in Figure 3. A design process proceeds through the Processes 0-8. However, most practitioners do not follow these separate steps in their daily practice, and very few of them can draw straightforward steps in any given process. Some researchers consider the process as an example of CPS (creative problem process), in which the entire process can be considered a spiral type. In the CPS domain, it places these processes on a circle and connects these processes with each other. Moreover, it indicates that all these processes are necessary, but the sequence can be varied. Based on these views and the experiences of consultants, [START_REF] Buijs | Modelling Product Innovation Processes, from Linear Logic to Circular Chaos[END_REF] proposes the circular version of the innovation model, as shown in Figure 4. All four stages of innovation (product launch and use, product development, strategy formulation and design brief formulation) are executed in the circular style. Leadership, culture, emotion, motivation, risk-taking and passion are the elements of the term "heart". Thus, we rely on the "heart," the competitive environment and the thinkingway the innovation team, as a way to start with out of the four stages.

Inspired by these analyses of Buijs, we conclude that the new version of the framework (Figure 3) similarly should be adapted as a circular style. Here, we see two reasons for the circular configuration in which we visualize our basic framework. Firstly, there is not any iteration in the basic framework.

Referring to the theories of innovation and the chained-linked model, we have already recognized innovation process is neither smooth nor linear, nor often well-behaved. Iterations and feedbacks are inevitable elements in the process. Secondly, due to the complexity andthe uncertainty of innovative design, it is impossible to finish the tasks we requested of each individual. So innovative design is by necessity compulsory teamwork, and different team members can perform these elementary processes at the same time. In other words, these processes are carried out in a parallel fashion. Therefore, our conception of a circular model with three layers provides an initial setting, into which we put these elementary processes of the basic framework. Figure 5 shows the circular process model. It consists of three layers: (1) Mapped onto in the middle layer are the nine elementary processes (described by Table2), with the iterations between these processes shown as dotted lines. The yellow parts denote that these processes are creative processes as in Figure 3, while the blue ones are not. The solid black lines represent the transformation and comparison between these design states (described by Table1).

(2) In the central layer, we transform the term "heart" into "designer", which drives these elementary processes.

(3) The outer layer focuses on the external environment that contains the influencing factors. Through the change of the form, the iterations and the parallel fashion can be well integrated. In the next section, the details regarding the centre and outer layers will be discussed.

Heart

Modifying the basic framework by introducing the designer and the environment space

In section 3.2, we configured the basic framework with a circular structure having three layers which shows the process of innovative design. However, there still are some questions: How does the designer drive these elementary processes? What is the external environment? How do the three layers interact with each other? In this section we propose answers to these questions. Design is a "reflective conversation with the materials of a design situation" [START_REF] Schön | The Reflective Practitioner: How Professionals Think In Action[END_REF]]. Through the perception of the design situation and the recollection of experiences, the designer decides on further actions (when to do what). That is to say, the designer should actively respond to their actions rather than simply perform the recommended actions. This means that the interaction between the designer and the external environment determines the design process. Furthermore, the considerations taken into account by the designer have already proven to be useful for improving innovation and design process. For example, it is useful to know which design activities are not able to have a positive effect in achieving the desired goal and which can improve the efficiency of design activities. Similarly, the designer plays a vital role in the development of innovation. Designers apply their skills, perspectives and domain knowledge to the innovation process based on their background and experience. Consequently, designers are generally recognized as an engine that drives the process, and in terms of the model a medium of interaction between different layers. The environment space is of significance in that it both influences and changes innovative design at particular moments. We define the environment space as a set of the external factors outside the designer. There are four reasons for introducing the environment space into innovative design. Firstly, innovative design should provide a means to partition the concepts and the knowledge in order to be more innovative [START_REF] Hatchuel | C-K design theory: an advanced formulation[END_REF]. It is impossible to acquire or generate new knowledge or concepts without the support of the environment. Secondly, the innovativeness of innovative design depends on the stakeholders involved. The stakeholders evaluate the value and the novelty of an innovative solution by communicating with the designer. According to the definition of the term "environment space", the stakeholders also should belong to one part of the environment. Thirdly, the process of innovative design is treated as a system, whichevolves by interacting with its environment. These environmental factors, such as organization architecture, organization cultureand organization strategy, influence the process of innovative design. Lastly, the designer is situated in the environment.The designer perceives and interprets design problems differently as environments change.

Then, a legitimate question that follows from this is how we consider the designer and the environment space together in the circular model. The point can be summarized by the term of reflection-in-action, which views that the entiredesign process as it proceeds from the designer"s perceptions. As we have discussed before, the term reflection-in-action is relevant to the notion of situatedness. [START_REF] Gero | The situated function-behaviour-structure framework[END_REF] describe situatedness as the interaction of the external world, the expected world and the interpreted world, see Figure 6. The external world, namely the environment space, is composed is represented as outside of the designer, while the interpreted word and the expected world are built up on the inside of the designer world. The interpreted world refers to sensory experiences, precepts and concepts of the designer. In the expected world, the actions that the designer will perform are imagined. These three worlds are dynamically coupled with each other through three types of processes: interpretation, focusing and action.

External World

Expected

World Interpreted World Focussing Action Interpretatinon

The internal world of Designer The environment space Figure 6. Situatedness as the interaction of three worlds Although Gero and Kannengiesser reconstruct these elementary processes of the FBS model based on the new framework, the external world is not clearly demarcated. That is to say, which elements does the external world include? According to the definition of the external world, the representations of these elements should be able to describe an external influence on the designer"s interpretation, and moreover, these are objective and independent of the designer. Furthermore, because the external world is dynamic with time, it should be seen as a time dependent variable. [START_REF] Reymen | A domain-independent descriptive design model and its application to structured reflection on design processes[END_REF] develop a concept called design situation to model design which defines the state of the combined concepts at a certain point in time. Here, we adopt the basic framework of the concept "design situation" to describe the external world, including the state of the design context, the state of the stakeholders, the design state and the state of the design process, as showed in Figure 7. The definitions of relevant concepts are listed in Table3. Having described situatedness and the concrete elements of the external world, we are now able to develop a more detailed framework of innovative design. In Figure 8, we show the three worlds of described situatedness, similar to Figure 6. Comparing with the initial circular model (Figure 5), the central layer "designer" is decomposed into the interpreted world and the expected world, and the middle layer and the outer layers consist of the external world. Through the pertinentdefinitions of Table 3, the design state and the state of the design process are located in the inner part of the external world, while the state of the design context and the state of the stakeholders are put into the outer part. In Figure 8, we represent more explicitly the process of innovative design as a dynamic environment, and also clearly depict the influencing elements.

Integrated descriptive model of innovative design

In the section we will attempt to map these elementary processes of the basic framework onto the situated model of innovative design. That is to say, it combines the Figure 5 and Figure 8 into an integrated descriptive model, as shown in Figure 9.

As can be seen, the layer shows the enhanced version of Vermaas and Dorst"s model based on the FBS model. The solid lines completely represent the different elementary processes (transformations and comparisons between these design states that are described in Table 1) and which can be found listed in Table 2. The Vermaas and Dorst" model focuses on the internals of the design process, but does not treat the innovative process, e.g. the process from requirement to purpose. Another change of the basic framework can be viewed as the adoption of the viewpoint of circularity in the basic framework. Because the complexity and uncertainty of innovative design, these elementary processes are impossible to be performed in sequence rather than in parallel fashion. Moreover, the feed forward loops and feedback loops incorporate the source of innovation. The circularity of these processes has been included in the structure in order to reflect these characteristics of innovation process.

Additionally, the integrated model consists of four layers, in which the and layers respectively correspond to the expected world and the interpreted world of the designer, and the and layers represent the external world of the designer. This is a consequence of introducing the designer and the environment space in order to deal with the designer"s interaction process with the external world and within oneself. In practice, the designer takes command of what he is "seeing" and "doing" in the design process, and communicates this information to the stakeholders, then makes actions in the current design context. The basic framework in Figure 3 is not able to reflect these interactive activities.

Action

Conclusion

In this paper, we first construct the basic framework of innovative design based on the FBS model, subsequently modify it by considering the innovation process, the designer and the environment space, and eventually propose a descriptive model of innovative design. The model reflects the interaction process between the designer, the design process and the outside environment.Whilst acknowledging that the process of innovative design observed in practice is more chaotic than the current representation suggests,understanding the relationships of the three parts can improving the process of innovative design. The advantage are summarised as follows: The model helps designers to better find the successful factors of innovative design, to integrate different aspect of an innovative design situation, to evaluate communication with the external environment. In addition, Insight into the layer model can reveal which processes may generate creative output. Firms can focus on these processes in order to enhance creative performance and the quality of the product designed. The current model is still a high-level description for innovative design.It can be a starting point for developing more effective management methods and tools to be used by firms in practices, especially product design. These methods and tools have the advantages that creativity activities are performed regularly based on a systematic approach. This should help to decrease the uncertainty and the complexity of innovative design, to generate more innovative ideas, and to design innovative products more quickly. Having a portfolio of detailed models and detailing available will help firms and designers to better deal with the process. Therefore, future studies should answer the question "How and to what level should we decompose the high-level model?" In addition, prescriptive studies should be conducted with designers to build up structured methods for different design processes.

Figure 2 .

 2 Figure 2. The chained-linked model

Figure 3 .

 3 Figure 3. The new basic framework of innovative design

  Figure4.The circular model of innovation Figure5.The initial circular model of innovative design

Figure 7 .

 7 Figure 7. The external world of the designer Figure 8. The situated model of innovative design

Figure 9 .

 9 Figure 9. The descriptive model of innovative design

Table 3 . The relevant definitions and reasons of the external world

 3 the accomplished design activities, and grounds them on previous knowledge. It can be seen as the external information of the designer.

	Concept	Definition	Why should it consider these concepts?
	Design	It is described by the sets of factors	An innovation process is an adaptive system which
	context	influencing the design state and the design process at a certain time.	evolves by reacting to its context. The contextual factors are one of the determinants of innovation.
	Stakeholder	A party that has an interest in	The stakeholders play a key role to evaluate the value
		innovative design.	and novelty of an innovative solution by communication
			with the designer.
	Design state	It refers to the key concepts of the FBS model, including Purpose	The four state variables are the temporal solutions of design. On the basis of the perceptions of these state
		Function, Behaviour and	variables, the designer makes expected actions.
		Structure(see Table 1)	Therefore, according to the definition of the external
			world, design states are one of parts of the environment.
	Design	It is defined as a series of design	The designer learns which activity were not successful

from