Marie Jardat

Vincent Dahirel

Fabrice Carnal

Dynamics of ions in model charged porous media: Influence of polyelectrolytes

We have studied by means of Brownian dynamics simulations the dynamics of small ions in model charged porous media. We have focused on the influence on this dynamics of short polyelectrolytes with the same charge sign as the solid phase. We have compared the selfdiffusion coefficients of counterions of four families of systems with different compositions (presence or not of charged obstacles and of polyelectrolytes). Our main result is that the presence of polyelectrolytes only modulates a behavior we had already observed before: In the domain where the porosity is relatively high compared to the Debye length, diffusion coefficients of counterions increase when the porosity decreases. Moreover, we have shown that, in the systems investigated here without charges on obstacles, the self-diffusion of counterions is mostly affected by the presence of polyelectrolytes and not by concentration effects, contrarily to co-ions. Also, we have seen that even if the density of the probability of presence of counterions in the vicinity of polyelectrolytes is much higher than that on obstacles, the residence time of ions around polyelectrolytes is shorter than around charged obstacles.

A C C E P T E D

M A N U S C R I P T

ACCEPTED MANUSCRIPT

I. INTRODUCTION

The dynamics of charged species in charged porous media is involved in various industrial and environmental situations, for instance in seawater desalinization processes or in the context of soil pollution by heavy metal ions. Numerical simulations have been used for many years to address this issue as they can bring microscopic explanations to observations or measurements at larger scales. For example, atomistic simulations of ions in clays are useful in the context of the confinement of toxic or radioactive waste in deep geological formations [1][2][3][4][5]. Simple models at intermediate scales can also be helpful. Partly-quenched systems, which consists of two components, a fixed one and an invading fluid, are representative of porous media imbibed by a fluid. In collaboration with V. Vlachy and B. Hribar-Lee, we have in the last decade applied replica Ornstein-Zernike calculations and Brownian dynamics simulations to understand the thermodynamic and transport properties of simple electrolyte solutions in partly-quenched systems. The solvent was in every case treated as a continuous medium. We have investigated the properties of monovalent and divalent ions in the presence of an electroneutral matrix of ionic obstacles representing disordered media [6]. We have also provided an explanation to the counterintuitive observed behavior of ions in some charged nanoporous media, where counterions of the porous medium are found to diffuse faster than co-ions (see for example [7,8]). In this case, we investigated the properties of monovalent ions in the presence of fixed macroions and found that the diffusion coefficient of counterions had a non-monotonic behavior as a function of the density of obstacles [9]. It increased slightly at small densities of obstacles, and had the expected decrease at high densities. Moreover, in dilute matrices, counterions were slower than co-ions, because they were strongly attracted by highly charged obstacles. In more concentrated matrices, where the average size of the pores was small, counterions were faster than co-ions because they could take an advantage of hopping from one obstacle to another, whereas the co-ions, repelled by obstacles, needed to diffuse in the crowded space between the layers of counterions surrounding the obstacles. The

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

influence of the electrostatic confinement of charged obstacles added to their steric confinement was also observed in minimal models of charged porous media where moving ions and obstacles had the same size [10,11].

Here, we would like to address the following new question. What happens if the invading solution also contains polyelectrolytes with a charge of same sign as the obstacles ? This generic situation may occur in soils in the presence of natural organic matter, which is known to play an important role in the mobility of cations (see for example [12][13][14][15]). It was shown recently that the polyacrylic acid (5100 g/mol)

represents a good analogue of humic acids in terms of size and charge in solution.

We recall that humic acids are components of the natural organic matter [START_REF] Roger | MANUSCRIPT loids and Surfaces A: Physicochem[END_REF]. At basic pHs, humic acids and PAA are negatively charged, so that in the presence of a negatively charged nanoporous medium, like clay, a competition between the attraction by the solid phase and by the polyelectrolytes would exist for cations.

Our aim here is thus to see if the presence of negatively charged polyelectrolytes modifies the main features of the dynamics of ions in nanoporous charged materials.

To answer this question, we propose to study simple partly-quenched systems by Brownian dynamics simulations.

The systems we are investigating consist of two components. The first one consists of large, negatively charged obstacles with a regular distribution. The second component consists of polyelectrolytes modelled as flexible chains of spherical monomers and of small ions. The whole system is assumed to be embedded in a continuous solvent modelling water. In what follows, we use Brownian dynamics simulations to compute the self-diffusion coefficients of ions in these systems. In this first study, we are not interested in the dynamic properties of the polyelectrolyte itself, even if polyelectrolytes are also mobile. More precisely, we investigate the influence of the distance between obstacles, which controls the porosity of the medium, on the ionic diffusion. In order to have a better insight into the main parameters which influence the dynamics of ions, we also run calculations (i) in the case where obstacles are uncharged, (i) in bulk solution with polyelectrolytes but without any obstacle, (iii) in the presence of charged obstacles but without any polyelectrolyte. We briefly

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

summarize the method and give the technical details in Section II, and we present the results in Section III.

II. MODEL AND METHODS

A. Brownian dynamics simulation

Brownian dynamics simulations based on an overdamped Langevin equation are used to compute the trajectories of moving ions. More precisely, the displacement of N particles included in the simulation box from time t to time t + δt reads [17]:

r(t + δt) = r(t) + βD • F(t)δt + R, (2.1)
where β = 1/k B T with T the temperature, D • is the self-diffusion coefficient of particles at infinite dilution, δt is the time increment, r is the 3N -dimensional configuration vector, and F is the total force acting on the particles at the beginning of the step. R is a random displacement, chosen from a Gaussian distribution with zero mean, R = 0, and variance RR T = 2D • δtI, where I is the identity matrix.

In simulations of charged particles the efficiency of the method is improved by evaluating the probability of each displacement using the smart Monte Carlo criterion [18,19].

Once the trajectories of particles are computed, the self-diffusion coefficient D of a given species is obtained from the autocorrelation function of the forces:

D = lim t→∞ D • - 1 3 t 0 (βD • ) 2 F(t 0 + t) • F(t 0 ) t 0 dt (2.2)
where F(t 0 ) and F(t 0 + t) are the forces on the particle at some arbitrary initial time t 0 and at some later time t 0 + t, respectively. Details on the simulation method applied to electrolyte solutions can be found in [19].

B. Model

The charged nanoporous medium is modelled as a cubic centered lattice of spheres which are much larger than moving ions. The simulation box thus contains two obstacles situated on the diagonal of the cubic box. It contains also five chains of 20 charged monomers each, representing the polyelectrolytes. The amount of small ions depends on the system, and there is in every case a number of monovalent counterions at least sufficient to neutralize the obstacles and the polyelectrolytes. In some cases, a +1 : -1 electrolyte is added in the simulation box, so that monovalent co-ions are also present. A typical snapshot of the simulation box is given in Fig.

II B.

The interaction pair potential between charged sites i and j (ions, obstacles or monomers) has a short-range and Coulomb part

U ij (r) = |Z i Z j |e 2 4πǫ 0 ǫ 1 12(a i + a j ) a i + a j r 12 + Z i Z j e 2 4πǫ 0 ǫ 1 r (2.3)
where ǫ 0 is the permittivity of the vacuum, ǫ is the relative permittivity of the solvent, e is the elementary charge, Z i is the valency and a i the size parameter of

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

particle i. r is the distance between centers of charge i and j. The choice of the shortrange part is such that the minimum of the interaction potential between particles of opposite charges is localized at the distance a i + a j . The model represents a "soft" version of the primitive model of electrolyte solutions, already used in related studies [6,9]. Polyelectrolytes are flexible bead-spring chains composed of uniformly charged monomers with the same physical size and magnitude of charge as salt ions.

There are 20 monomers per chain. The bonding potential between neighboring beads along the polymer backbone, V FENE (r), is given by the finitely extensible nonlinear elastic (FENE) potential,

V FENE (r) = - 1 2 k FENE R 2 0 ln 1 - r 2 R 2 0 (2.4)
where k FENE is the spring constant and R 0 is the maximum extension of the bond.

In this study, the obstacles have a size parameter a obstacle = 2.0 nm with a charge Z obstacle = -50, the ions of the electrolyte (counterions of obstacles and of polyelectrolytes and ions of the added salt) have a size parameter a ion = 0.15 nm with a charge Z ion = ±1, monomers of the polyelectrolyte have a size parameter a mono = 0.15 nm with a charge Z mono = -1. The parameters of the bonding potential are k FENE = 30k B T /(2a mono ) 2 and R 0 = 3a mono .

We study four families of systems. The first one, described above, contains charged obstacles, polyelectrolytes and salt ions; We refer to this family of systems as OB-PE (for OBstacles-PolyElectrolyes). In the second family, obstacles are neutral;

We refer to this family as neutralOB-PE. The third family of systems corresponds to a bulk solution of polyelectrolytes without any obstacle; We refer to this family as BulkPE. In the last one, the electrolyte solution is embedded in charged obstacles but there are no polyelectrolytes; We refer to this family as OB.

In what follows, we use the effective polyelectrolyte concentration in the presence of obstacles c eff , defined by c eff = N PE /V free where N PE is the number of polyelectrolytes in the simulation box, and V free is the free volume, i.e. the volume of the simulation box decreased by that occupied of obstacles (V free = L box 3 -8 3 πa 3 obstacle ).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

We will also refer to the Debye length λ D , which we calculate as follows:

λ D = ǫ 0 ǫk B T V free N ion e 2 1/2 (2.5)
with N ion the total number of small ions in the simulation box, which means that we neglect the contribution of charged obstacles and of polyelectrolytes to the screening length.

In the OB-PE systems, 2 obstacles and 5 polyelectrolytes are put in the simulation box, so that there are 100 counterions. For each family of systems, the size of the simulation box takes four different values, which allows us to vary the porosity of the system. In systems OB-PE, neutralOB-PE and OB, the box length varies between L box = 15.0 nm, and L box = 7.90 nm. The distance d obst between the solid surfaces of obstacles varies thus between 9 nm for the largest and 2.8 nm for the shortest. In every case, we are thus in regimes where a polyelectrolyte chain, whose maximum extension is about 9 nm and whose bead diameter is 0. Note that the size of the symbols in Figs. 4 is of the order of the uncertainty.

III. RESULTS AND DISCUSSION

We show in Fig. 2 the self-diffusion coefficients of counterions (plain line) and of co-ions (dashed line) in systems which contain neutral obstacles and polyelectrolytes for the largest distance between obstacles, i.e. the largest porosity (system neutralOB-PE1), and for the smallest one (system neutralOB-PE4). More precisely effective concentration of polyelectrolytes varies by a factor 6 between these two systems but it remains low in both cases (less than 0.02 mol/L).

It is striking in Fig. 2 to see that the diffusion of co-ions is much more influenced by the variation of the porosity than that of counterions in this kind of systems.

The diffusion coefficient of co-ions is decreased by about 10 % when the porosity decreases, whereas that of counterions stays almost the same. As for co-ions, they are repelled by polyelectrolytes and move in the free space between obstacles and polyelectrolytes. The number of small ions in the simulation box is the same in both systems, so that their concentration varies when the porosity vary. They are

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

actually ten times more concentrated than polyelectrolyte chains in the free space (about 0.2 mol/L in neutralOB-PE4). The variation of their concentration between systems neutralOB-PE1 and neutralOB-PE4 by a factor 6 explains the decrease of their self-diffusion coefficient with the porosity.

As for counterions, the situation is different: They are attracted by the polyelectrolytes, and in average, a constant amount of them stays in the vicinity of the chains. This ionic condensation of counterions on the polyelectrolyte chains is illustrated in Fig. 3 which shows the radial distribution functions between counterions and PE beads. The rdfs present a high peak which decreases when the whole system is more concentrated (compare dashed and plain lines, respectively). We also see in Fig. 3 that the g(r) is almost not affected by the presence of neutral obstacles.

Finally, we see in Fig. 2 that the diffusion of counterions tends to be slowed down when the concentration increases, but this effect is very small. Indeed, as the polyelectrolytes diffuse about 20 times slower than small ions, because they consists of 20 beads, the mean diffusion coefficients of counterions is strongly decreased compared to the infinite dilution. But, as the concentration of the polymeric chains is relatively low in every case, their diffusion is not considerably affected by concentration effects. In conclusion, we can say that in neutralOB-PE systems, the self-diffusion of counterions is mostly affected by the presence of polyelectrolytes and not by concentration effects, contrarily to co-ions. In what follows we focus on the diffusion coefficients of counterions.

We give in Fig. 4 the values of the ratios D/D • for the three families of systems which contain obstacles, namely OB-PE, neutral-OB-PE and OB, as a function of the distance between the surfaces of obstacles. We recall that at a given distance between obstacles, the Debye length is the same in the three different systems. First of all, we recover for systems neutralOB-PE the results already presented in Fig. (systems OB, black triangles) follow the trends already observed in our previous studies [9]. In this regime of relatively low concentration of obstacles, the selfdiffusion coefficient of counterions increases when the distance between obstacles decreases, i.e. when the porosity decreases, from right to left in Fig. 5. We had interpreted this counterintuitive behavior in previous studies as an effect of electrostatic attractions between charged obstacles and counterions. When the distance between obstacles decreases, the Debye length also decreases, so that the electrostatic condensation of counterions on obstacles is less intense. This trend is kept in the presence of polyelectrolytes (systems OB-PE, red circles): In this case also we observe an increase of the diffusion coefficient of counterions when the porosity decreases (from right to left). The radial distribution functions (rdfs) between counterions and charged obstacles are presented in Fig. 5: They confirm that counterions are less condensed on charged obstacles when the distance between obstacles decreases (compare dashed and plain lines in both cases with and without polyelectrolytes), as the peak of the rdf decreases.

The comparison of black and blue curves in Fig. 4 shows that the diffusion of counterions is faster when they are attracted only by obstacles (systems OB) than when they are attracted only by polyelectrolytes (systems neutralOB-PE), even if polyelectrolytes are mobile, contrarily to obstacles. This is due to an attraction of counterions by polyelectrolytes stronger than by obstacles (see the g(r) in Fig. and 5), which implies a smaller amount of free (and faster) counterions. Finally, if we compare red and blue plots in Fig. 4, we see that for short distances between obstacles, counterions are faster when the obstacles are charged. But, for large distances between obstacles, counterions are faster when the obstacles are neutral.

The radial distribution functions displayed in Fig. 5 also show that at a given porosity, the presence of polyelectrolytes in the system decreases the peak between counterions and charged obstacles (compare both dashed lines for example). In the presence of charged polymers, there is indeed a competition between the attraction by charged obstacles and by polyelectrolytes. This illustrates well the competition between polyelectrolytes and obstacles effects. time of course depends on the definition of the "vicinity", which corresponds actually to a cut-off distance under which we assume that counterions are condensed on the bigger charged object. We show in Fig. 6 the evolution of the mean first exit time as a function of the cutoff distance in system OB-PE2. The first exit time is computed as follows. Each counterion in the simulation box can belong only to one condensed layer of one larger charged object (if it belongs to several ones, we kept the layer which corresponds to the shortest distance). We record for each trajectory the time a counterion spends in a given layer from the beginning of the trajectory until its exit. This quantity is averaged over all initially condensed counterions. As it can be seen in Fig. 6, for short cutoff distances, this time is short, because the distance to be covered to exit the layer is short. The mean first exit time first increases with the cutoff distance and then decreases, because for long cutoff distances, the number of counterions only slightly attracted to the charged object and initially situated close to the cutoff distance increases. These counterions have short exit times, being close to the boundary. We have plotted the mean first exit time as a function of the cutoff distance for the three families of systems OB-PE, neutralOB-PE and OB. We have then kept the mean first exit time which corresponds to the first peak of this function. Results are collected in Fig. 7 for the residence time on polyelectrolytes, and Fig. 8 for the residence time on charged obstacles.

In every case, the exit time increases with the distance between obstacles, which means that counterions spend more time in the vicinity of polyelectrolytes or of charged obstacles when the whole concentration of the system decreases. The stronger the ionic condensation, the longer the time spent in the condensed layer.

As for the condensation around polyelectrolytes, we obtain almost the same values with neutral and with charged obstacles, which are rather small. The exit time is slightly larger when obstacles are uncharged, which is expected because in this case there is no competition for the attraction of counterions. Nevertheless, we showed above that, in the presence of neutral obstacles, no significant variation of selfdiffusion coefficients with the distance between obstacles appeared. This can be explained by the fact that the exit time is short in every case. In other words, the number of counterions condensed on polyelectrolyte chains is large in average but with a fast dynamic of exchange. The mean first exit time on charged obstacles displayed in Fig. 7 is larger than that on polyelectrolytes, and it decreases when there is a competition due to the presence of polyelectrolytes (black plot compared to the red one). Finally, we see here that, even if the average condensation on charged

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

obstacles is lower than on polyelectrolytes, the residence time is larger, which means that there are less exchanges within the condensed layer.

IV. CONCLUSION

Our aim here was to see if the presence of negatively charged polyelectrolytes modifies the main features of the dynamics of ions in model nanoporous charged materials. In this case, counterions are not only attracted by the solid phase (here modelled as negatively charged spheres on a regular lattice) but also by the negatively charged electrolytes. We have compared the self-diffusion coefficients of counterions in four families of systems which differ from their composition (presence or not of charged obstacles and of polyelectrolytes). Our main result is that the presence of counterions does not change the main features of what we observed before:

Diffusion coefficients of counterions increase when the distance between charged obstacles decreases, even in the presence of a low concentration of polyelectrolytes.

Moreover, we have shown that, in the systems investigated here without charges on obstacles, the self-diffusion of counterions is mostly affected by the presence of polyelectrolytes and not by concentration effects, contrarily to co-ions. Also, we have seen that even if the density of the probability of presence of counterions in the vicinity of polyelectrolytes is much higher than that on obstacles, the residence time of ions around polyelectrolytes is shorter than around charged obstacles.

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Highlights

• The diffusion of ions in model charged porous media is studied.

• Counterions diffuse faster when the porosity decreases.

• Counterion diffusion is mostly affected by polyelectrolytes.

• Ions spend less time around polyelectrolytes than around obstacles.

FIG. 1 :

 1 FIG. 1: Snapshot of a typical simulation box: obstacles are in blue, polyelectrolytes in red and small ions in grey. Obstacles are totally included in the cubic box so that they are actually distributed on the nodes of a cubic centered lattice with a lattice constant a = 2 L box .

  3 nm should be able to diffuse in between obstacles. The BulkPE systems have slightly different box lengths in order that the polyelectrolyte concentration is exactly the same as the effective one in the presence of obstacles. Moreover, in systems without obstacles, in those without polyelectrolytes and in those with neutral obstacles, we have added a salt such that the Debye length keeps the same value as in the presence of obstacles. More precisely, we have thus added 50 cations (Z = +1) and 50 anions (Z = -1) in families of systems neutralOB-PE, OB and BulkPE, which corresponds to a contribution to the Debye length equivalent to that of the 100 counterions in family OB-PE. The Debye length is in every case shorter than the distance between the surfaces of obstacles. The systems are labelled with indices 1 to 4 in each family, with the index taking increasing values as the box length decreases. The main characteristics of the studied systems are summarized in Table I. The self-diffusion coefficients of ions are averaged over 10 to 20 independent runs of total duration greater than the typical time scale needed to diffuse over the whole simulation box. In what follows, the results are expressed as the ratios D/D • between the computed self-diffusion coefficients and the value D • at infinite dilution. A C C E P T E D M A N U S C R I P T

FIG. 2 :

 2 FIG. 2: Self-diffusion coefficients of small ions divided by the value at infinite dilution as a function of time in the presence of neutral obstacles and of polyelectrolytes (systems neutralOB-PE). Counterions (Z ion = +1): plain line; Co-ions (Z ion = -1): dashed lines. The values obtained in system neutralOB-PE1 (blue) correspond to the largest distance between the surfaces of obstacles (d obst ≃ 9 nm). The values obtained in system neutralOB-PE4 (red) correspond to the smallest distance between the surfaces of obstacles (d obst ≃ 3 nm)

2 :FIG. 3 :

 23 FIG. 3: Radial distribution functions between counterions and polyelectrolytes. Systems neutralOB-PE1 and neutralOB-PE4: black lines, respectively plain and dashed lines. Systems BulkPE1 and BulkPE4: red lines, respectively plain and dashed lines.

FIG. 4 :

 4 FIG. 4: Self-diffusion coefficients of counterions as a function of the distance between the surfaces of obstacles for the three different families of systems. Red filled circles: Systems OB-PE, which contain charged obstacles and polyelectrolytes. Black filled triangles: Systems OB, which contain only charged obstacles. Blue filled diamonds: Systems neutralOB-PE, which contain neutral obstacles and polyelectrolytes.

3

 3

FIG. 5 :

 5 FIG. 5: Radial distribution functions between counterions and charged obstacles. SystemsOB-PE1 and OB-PE4: red lines, respectively plain and dashed lines. Systems OB1 and OB4: black lines, respectively plain and dashed lines.

FIG. 6 :

 6 FIG. 6: Mean first exit time of counterions from the vicinity of polyelectrolytes (in black)and from the vicinity of charged obstacles (in red) for system OB-PE2 as a function of the cutoff distance under which counterions are assumed to be in the condensed layer (or in the "vicinity").

FIG. 7 :

 7 FIG. 7: Mean first exit time of counterions from the condensed layer around polyelectrolytes as a function of the distance between obstacles. Systems OB-PE: red circles; Systems neutralOB-PE: black diamonds.

FIG. 8 :

 8 FIG. 8: Mean first exit time of counterions from the condensed layer around charged obstacles as a function of the distance between obstacles. Systems OB-PE: red circles; Systems OB: black triangles.

TABLE I :

 I NameL box /nm d obst /nm λ D / nm c eff / mmol.L -1 Main characteristics of the studied systems.

	ACCEPTED MANUSCRIPT
	OB-PE1	15.0	8.99	1.35	2.5
	OB-PE2	11.0	5.53	0.84	6.6
	OB-PE3	9.7	4.40	0.68	9.9
	OB-PE4	7.9	2.84	0.49	19.6
	neutralOB-PE1	15.0	8.99	1.35	2.5
	neutralOB-PE2	11.0	5.53	0.84	6.6
	neutralOB-PE3	9.7	4.40	0.68	9.9
	neutralOB-PE4	7.9	2.84	0.49	19.6
	BulkPE1	14.9		1.35	2.5
	BulkPE2	10.8		0.84	6.6
	BulkPE3	9.5		0.68	9.9
	BulkPE4	7.5		0.49	19.6
	OB1	15.0	8.99	1.35	0
	OB2	11.0	5.53	0.84	0
	OB3	9.7	4.40	0.68	0
	OB4	7.9	2.84	0.49	0

Acknowledgements: The authors thank the financial support from the French National Agency for Research (ANR) under grant ANR-09-JCJC-0082-01. MJ also thanks Vojko Vlachy for the interesting and enjoyable time they spent together in Paris and in Ljubljana discussing not only the results they obtained, but also what