Antoine Carof

Virginie Marry

Mathieu Salanne

Jean-Pierre Hansen

Pierre Turq

Benjamin Rotenberg
email: benjamin.rotenberg@upmc.fr

Coarse graining the dynamics of nano-confined solutes: the case of ions in clays

Keywords: Generalized Langevin Dynamics, confinement, memory function, clays, cesium

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The dynamics of confined fluids differs dramatically from their bulk counterparts, especially when the length scale of the confinement becomes comparable to the molecular size [START_REF] Bocquet | Hydrodynamic properties of confined fluids[END_REF]. In such cases, the hydrodynamic flows may depart from the predictions of continuous theories, even when modified boundary conditions at the wall are introduced, because of the layering due to the finite size of the fluid molecules [START_REF] Bocquet | Flow boundary conditions from nano-to micro-scales[END_REF]. Under even more extreme confinement, where one or two molecular layers of fluid are present, hydrodynamic descriptions are not relevant anymore. While molecular simulation remains the reference theoretical tool to investigate such situations, resort to continuous models is obviously necessary if one wishes to understand the evolution of large systems over long time scales. For example, diffusion models may still be relevant, since the motion of fluid molecules on very long time scales remains diffusive. Nevertheless such models may not accurately reflect the underlying microscopic mechanism, which may differ significantly from that for bulk diffusion. In particular, while the Brownian theory of diffusion, which relies on a separation of time scales between the solute and its surrounding medium is expected to break down [START_REF] Bocquet | High friction limit of the Kramers equation: The multiple time-scale approach[END_REF][START_REF] Bocquet | From a stochastic to a microscopic approach to Brownian motion[END_REF].

As an example of extremely confined solute, we investigate the case of cesium ions confined in the interlayer porosity of montmorillonite clays, with a very low water content (a single water layer with 6 water molecules per ion), for which previous attempts to build a continuous solvent model failed to reproduce the dynamics predicted by molecular simulations [START_REF] Rotenberg | A multiscale approach to ion diffusion in clays: Building a two-state diffusion-reaction scheme from microscopic dynamics[END_REF][START_REF] Dufrêche | Bridging molecular and continuous descriptions: the case of dynamics in clays[END_REF]. We first derive a Generalized Langevin Equation (GLE) for the motion of a solute in an external potential using the Mori-Zwanzig formalism. We then compute the corresponding memory function from molecular simulations. The purpose of the present work is not to perform GLE simulations using this memory kernel, but rather to investigate its properties. In particular, we discuss the separation of time scales between the confined ions and solvent. We show that such a separation is not achieved and analyze the relative contribution of the ion-surface, ion-solvent and ionion interactions to the dynamics. This explains why the dynamics of this ion under extreme confinement cannot be described by a simple Langevin equation, as was previously attempted.

2 Theory: from Newtonian to Langevin and Generalized Langevin dynamics

Continuous solvent model for ions confined in clay interlayers

We consider here the case of ions confined in the interlayer porosity of clay minerals. These layered aluminosilicate minerals are abundant in the Earth's crust and play an important role in many environmental and industrial processes. Their low porosity and their surface charge control the mobility of ions through their multi-scale porosities and the sorption of these ions on the various types of surfaces they offer. This in turn explains their consideration e.g. for the retention of toxic and radioactive waste. Experimental evidence suggests that in clay rocks the major pathways for cations is the so-called interlayer porosity [START_REF] Glaus | Diffusion of 22 Na and 85 Sr in Montmorillonite: Evidence of Interlayer Diffusion Being the Dominant Pathway at High Compaction[END_REF], where the counter-ions compensate the negative charged of the mineral layers (see Figure 1). The interlayer space may also contain some water, in quantities which depend among other on the surface charge density and the nature of the counter-ion. In particular, for montmorillonite clays with cesium counterions at low relative humidity (and/or high confining pressure), the most stable state corresponds to a water monolayer, with only ≈ 6 water molecules per cation. The distance between the confining surfaces is then comparable to the molecular and ionic size. The motion of ions confined in these sub-nanometer pores has been studied extensively by molecular simulation [START_REF] Marry | Microscopic simulation of structure and dynamics of water and counterions in a monohydrated montmorillonite[END_REF][START_REF] Marry | Microscopic Simulations of Interlayer Structure and Dynamics in Bihydrated Heteroionic Montmorillonites[END_REF][START_REF] Malikova | Na/Cs montmorillonite: temperature activation of diffusion by simulation[END_REF][START_REF] Malikova | Diffusion of Water in Clays on the Microscopic Scale: Modeling and Experiment[END_REF]. It is almost two-dimensional and diffusive in the long time limit, with a linear growth of the mean-square displacement [START_REF] Marry | Water diffusion in a synthetic hectorite by neutron scattering beyond the isotropic translational model[END_REF]. The corresponding diffusion coefficient is ≈ 10 times smaller than in a bulk, infinitely diluted aqueous solution. In a previous work, we have investigated the possibility to describe this confined motion by a continuous solvent model in which the ion evolves in the Potential of Mean Force (PMF), which captures the average effect of all its environment (clay surfaces, water molecules and other ions). The effect of a thermalizing bath is then accounted for via a constant and uniform friction and random forces [START_REF] Rotenberg | A multiscale approach to ion diffusion in clays: Building a two-state diffusion-reaction scheme from microscopic dynamics[END_REF][START_REF] Dufrêche | Bridging molecular and continuous descriptions: the case of dynamics in clays[END_REF]. The resulting model corresponds to simple Langevin Dynamics (for simplicity we write it here in the one-dimensional case):

dv dt = F P M F (t) m -γv(t) + R(t) m , (1)
where v is the velocity of the ion, m its mass, F P M F the mean force acting on the ion at its current position, γ the friction and R a random force. The distribution of the latter is taken as Gaussian with a zero mean and a variance R(t)R(t ) = 2γmk B T δ(t -t ), k B being the Boltzmann constant, T the temperature and δ the Dirac distribution. In general, Lanvegin dynamics refers to numerical simulations based on the coupled stochastic equations for the evolution of N ions interacting via the instantaneous force due to the surfaces and to the other ions. In the following, however, we will use Langevin dynamics to refer to the effective one-body dynamics described by Eq. (1). Figure 2 compares the Velocity Autocorrelation Function (VACF) from Molecular Dynamics (MD) simulations with the results of Langevin dynamics in the PMF extracted from molecular simulations, for several values of the friction coefficient. It shows that no choice of friction allows to capture simultaneously the initial decay, the depth and the position of the minimum [START_REF] Rotenberg | A multiscale approach to ion diffusion in clays: Building a two-state diffusion-reaction scheme from microscopic dynamics[END_REF][START_REF] Dufrêche | Bridging molecular and continuous descriptions: the case of dynamics in clays[END_REF].

0 1 2 3 t ( ps ) -1 0 1 2 C(t) ( Å 2 / ps ) MD γ = 5 ps -1 γ = 2 ps -1 γ = 1 ps -1
Figure 2. Velocity Autocorrelation Function (VACF) obtained by Molecular Dynamics (black curve) and by Langevin Dynamics using the potential of mean force extracted from molecular simulations, for different values of the friction γ (red, green and blue curves). None of these frictions is able to capture simultaneously the initial decay, the depth and the position of the minimum. Figure adapted from Ref. [START_REF] Dufrêche | Bridging molecular and continuous descriptions: the case of dynamics in clays[END_REF].

An underlying assumption of the Langevin equation ( 1) is the decoupling of time scales between the solute (ion) and its environment, which may not be achieved under extreme confinement conditions, despite the larger mass of the cesium ion compared to that of water molecules. The traditional way beyond this approximation is to describe the motion of the solute via a Generalized Langevin Equation (GLE), which in the absence of an external potential reads [START_REF] Hansen | Theory of Simple Liquids[END_REF]:

dv dt = - t 0 K(u)v(t -u) du + R(t) m . (2)
This equation introduces a memory kernel linked to the random force (noise) via the fluctuationdissipation relation:

K(t) = R(0)R(t) mk B T (3)
Note that the simple Langevin equation corresponds to an infinitely fast relaxation of the noise, compared to the evolution of the velocity, i.e. K(t) = 2γδ(t), with the friction coefficient:

γ = ∞ 0 K(t)dt. (4)
In the limit of an infinitely heavy solute, this friction coefficient can be computed as the integral of the force autocorrelation function (FACF) instead of the memory kernel K [START_REF] Mazur | Molecular theory of Brownian motion[END_REF][START_REF] Espanol | Force autocorrelation function in Brownian motion theory[END_REF]. In the following, we discuss the extension of this framework to include the interactions with the confining surfaces. We first derive the relevant quantities to be determined using the Mori-Zwanzig formalism of projection operators. We then compute these quantities from molecular dynamics simulations. We finally discuss the applicability of such a strategy in the present case.

Derivation of a Generalized Langevin Equation

The Mori-Zwanzig formalism provides a natural route from the molecular Newtonian dynamics to coarse-grained models in which the dynamics is projected onto a reduced number of variables.

We derive here a GLE for the evolution of the velocity of a solute under extreme confinement.

In order to set the stage, let us first recall the standard derivation of Eq. ( 2) in the absence of such a potential. The starting point is the Liouville equation for the velocity of the solute:

dv dt = iLv (5)
where the Liouvillian (super)operator L is related to the total Hamiltonian H of the system by: iL = {H, .}. Since we aim at reducing the description of the whole system to the sole velocity v of an ion, we introduce the two operators P and Q defined by their action on any observable property A as:

P A(t) = v(0)A(t) v(0) 2 v(0) = m v(0)A(t) k B T v(0) QA(t) = A(t) -P A(t) . (6)
The effect of P is to "project" the time-dependence of A on v, according to its correlation with the latter. Formal integration of the Liouville equation ( 5) allows to express the rate of change in the velocity from a given initial condition as:

dv dt = e iLt iLv(0) = e iLt F (0) m + e iLt Q F (0) m , (7)
where we have introduced the fact that P + Q = 1 and F (0) = m × iLv(0) the force acting on the solute arising from its interactions with the solvent. Using the parity properties of correlation functions, the hermitian properties of the Liouville operator and Dyson's relation [START_REF] Zwanzig | Nonequilibrium Statistical Mechanics[END_REF], one obtains the GLE (2) with the memory kernel (3) and the explicit formal expression of the noise:

R(t) = e iQLt F (0) . (8)
We now turn to the more complex case where an external potential acts on the solute, and separate the total force as F tot = F ext + F bath . Note that such a decomposition is very general and the choice of how to separate the force due to the external potential F ext and to the thermal bath F bath will be discussed below. Following the same procedure, we obtain the GLE:

dv dt = F ext (t) m - t 0 v(t -u)K(u) du + R(t) m (9)
Despite a very strong ressemblance with Eq. ( 2), this result does not simply correspond to adding the external force to the equation of motion. Indeed, the memory K and noise R are now obtained by considering the action of the projection operators on the bath force only, instead of the total force, with the final result:

K(t) = F tot (0)R(t) mk B T (10)
R(t) = e iQLt F bath (0) . (11)
Note that the external potential is also present in the memory function K. Of course, in the absence of an external potential this reduces to the previous expressions.

Comparing the exact result of Eqs. [START_REF] Marry | Microscopic Simulations of Interlayer Structure and Dynamics in Bihydrated Heteroionic Montmorillonites[END_REF][START_REF] Malikova | Na/Cs montmorillonite: temperature activation of diffusion by simulation[END_REF][START_REF] Malikova | Diffusion of Water in Clays on the Microscopic Scale: Modeling and Experiment[END_REF] to the model of Ref. [START_REF] Dufrêche | Bridging molecular and continuous descriptions: the case of dynamics in clays[END_REF] described in the previous section, we now can explicit the approximations underlying the latter. On the one hand, it identifies the external force with the Potential of Mean Force F P M F = -∇V P M F . This seems reasonable since the interaction with the clay indeed corresponds to an external potential, however the PMF also includes the average effect, for a given position of the solute, of the solvent and of the other solutes. On the other hand, it assumes that the remainder of the force F bath = F tot -F P M F relaxes very rapidly compared to the velocity of the solute so that its effect can be captured by a memory-less friction and a corresponding white noise. This assumption can be challenged by computing the memory kernel K from molecular simulations and comparing its relaxation with that of the VACF, as we now proceed to show. In addition, we analyze the effect of confinement on the separation of time scales between the solute and its environment by comparing the case of the cesium ion in the interlayer of montmorillonite and in bulk water. Finally, we discuss the contributions of the ion-clay, ion-water and ion-ion interactions to the total force and its relaxation.

3 Molecular simulations: methods

System

The simulated clay system, illustrated in Figure 1 The water molecules are described according to the SPC/E force field [START_REF] Berendsen | The missing term in effective pair potentials[END_REF], the cesium ions by the force field from Koneshan et al. [START_REF] Koneshan | Solvent Structure, Dynamics, and Ion Mobility in Aqueous Solutions at 25 C[END_REF], and the clayFF force field is used for the clay layers [START_REF] Cygan | Molecular Models of Hydroxide, Oxyhydroxide, and Clay Phases and the Development of a General Force Field[END_REF]. The latter are considered rigid and the corresponding atomic positions are set to the experimental data [START_REF] Maegdefrau | Die Kristallstruktur des Montmorillonits[END_REF]. Long-range electrostatic interactions are computed using Ewald summation [START_REF] Frenkel | Understanding Molecular Simulations, From Algorithms to Applications[END_REF]. Molecular dynamics in the NVT-ensemble are performed using a Nosé-Hoover thermostat with a time constant of 0.5 ps. The system is first randomized at 1000 K for 10 ps, then equilibrated at 298 K for 80 ps. The properties are subsequently determined in NVE runs with different lengths: 10 simulations of 35 ps and 4 simulations of 5 ns for the short and long time dynamics, respectively. Newton's equation of motion are solved using the Leap-Frog Verlet algorithm with a time-step of 1 fs. All simulations were performed using the DL POLY simulation package [START_REF] Smith | DLPOLY2 user manual Daresbury Laboratory[END_REF].

Potential of Mean Force

The two-dimensional PMF is computed from the equilibrium one-body density as: V P M F (x, y) = -k B T ln ρ(x, y). The latter is determined by averaging over ten trajectories of 35 ps, over the 48 ions and exploiting the hexagonal symmetry of the cavities on the surface. Since the mean force needs to be evaluated at any position in order to extract the bath force, we further use an analytical approximation to represent it using the following expression, which corresponds to the above-mentioned symmetry of the unit-cell:

V analy P M F (x, y) = -k B T 6 nx=0 6 ny=0 A(n x , n y ) cos 4πn x x a x cos 4πn y y a y + cos 2πn x (x + y √ 3) a x cos 2πn y (y -x √ 3) a y + cos 2πn x (x -y √ 3) a x cos 2πn y (y + x √ 3) a y , (12)
with a x and a y the dimensions of the unit cell. The coefficients A(n x , n y ) are determined numerically so as provide the best fit to the molecular simulation results, using the MINUIT software [START_REF] James | MINUIT -System for function minimization and analysis of parameter errors and correlations[END_REF]. The mean force F P M F is then computed at each time step as the derivative of this expression, evaluated at the position of the cesium ion.

Computation of the noise and of the memory kernel

While the extraction of the VACF and FACF from a MD trajectory is straightforward, as the velocity and force acting on each particles are readily available, the determination of the noise R, hence of the memory kernel K, is much more involved. Indeed, the projected forces do not evolve simply according to the Liouvillian as in Eq. ( 5), but according to the projected operator iQL. By differentiating the definition Eq. ( 11) with respect to time and using the properties of the projector Q, we obtain the following evolution equation for the noise [START_REF] Darve | Computing generalized Langevin equations and generalized Fokker-Planck equations[END_REF]:

dR dt = iLR(t) + F tot (0)R(t) k B T v (13)
which has to be solved for the initial condition:

R(0) = F bath (0) . (14)
The numerical algorithm allowing to reconstruct a posteriori the noise from a MD trajectory will be presented in Ref. [25], which discusses two algorithms based on different points of view (looking forward or backward in time) and illustrates their relative merits on the case of a Lennard-Jones fluid. For each configuration of the trajectory, the bath force is determined by substracting the analytical representation of the mean force from the total force and used as the initial condition to integrate Eq. ( 13). The memory kernel is finally obtained from its definition Eq. (10).

Force Splitting

In order to assess the relative contributions to the overall force on the confined ion due to the ionclay (C), ion-water (W) and ion-ion (I) interactions, these forces must be computed separately at each time step of the trajectory. The timescales on which they evolve and their possible coupling can then be quantified by their auto-and cross-correlation functions.

4 Results and Discussion

Potential of Mean Force

The two-dimensional PMF for a Cs + ion in a monohydrated montmorillonite is reported in Figure 3. It displays a global minimum corresponding to the ion above the center of the hexagonal cavity as well as secondary minima above the surface silicon atoms (i.e. coordinated by three surface oxygen atoms to which the silicon is bonded). Both situations can be observed on Figure 1. The diffusion of the Cs + ion from one cavity to another involves hopping through the secondary minima, with an activation barrier of ∼ 3k B T . These results are in qualitative agreement with those of Ref. [START_REF] Rotenberg | A multiscale approach to ion diffusion in clays: Building a two-state diffusion-reaction scheme from microscopic dynamics[END_REF][START_REF] Dufrêche | Bridging molecular and continuous descriptions: the case of dynamics in clays[END_REF] obtained with a different force field.

Short time dynamics

The main purpose of the present paper is to assess the relevance of a continuous solvent model to describe the dynamics of ions confined down to the molecular scale. The VACF for Cs + confined in the interlayer of montmorillonite in the monolayer hydration state is reported in Figure 4.

In agreement with previous simulations of the same system using a different force field [START_REF] Rotenberg | A multiscale approach to ion diffusion in clays: Building a two-state diffusion-reaction scheme from microscopic dynamics[END_REF], we find that the VACF decreases from k B T /m to a negative minimum value (≈ -0.25k B T /m) in approximately 0.5 ps, before decaying to 0 for t ≈ 1 ps). As mentioned above, the simple Langevin equation ( 1) has been shown to fail to reproduce this VACF obtained by molecular simulation [START_REF] Rotenberg | A multiscale approach to ion diffusion in clays: Building a two-state diffusion-reaction scheme from microscopic dynamics[END_REF][START_REF] Dufrêche | Bridging molecular and continuous descriptions: the case of dynamics in clays[END_REF]. We now investigate in more detail the validity of the main approximation underlying this model, namely that the thermal bath which models the deviation of the instantaneous force from the mean force, evolves on a faster time scale than the solute. In the limit of an infinitely heavy solute (compared to the solvent) the friction coefficient of the Langevin equation ( 1) is given by the integral of the force autocorrelation function (FACF) [START_REF] Mazur | Molecular theory of Brownian motion[END_REF][START_REF] Espanol | Force autocorrelation function in Brownian motion theory[END_REF]. Such an assumption may seem reasonable in the case of the Cs + which is ≈10 times heavier than a water molecule. The FACF, also reported in Figure 4 decays to zero more rapidly than the VACF, within 0.5 ps, after a significantly faster initial decrease (0.1 ps) to a negative minimum (≈ -0.2 F 2 ) followed by a shallower one (≈ -0.1 F 2 at 0.3 ps). Under the above approximation, the relatively fast decorrelation of the force compared to the velocity, even though not complete, may suggest that the Brownian approximation is reasonable.

However, the relevant correlation function to define the memory kernel K is the noise autocorrelation function (NACF), as can be seen in Eq. [START_REF] Malikova | Na/Cs montmorillonite: temperature activation of diffusion by simulation[END_REF]. The NACF, computed from the noise extracted from the trajectory as explained in Section 3.3, is reported for the confined ion in Figure 4. It coincides with the FACF at very short times (up to 50 fs) and its variations up to 250 fs, with a miminum followed by a maximum, clearly ressemble that of the FACF in this intermediate regime. Nevertheless, there are significant differences between the memory kernel and the FACF. Firstly, contrary to the latter, the memory is always positive. Secondly, and most strinkingly, the memory appears to level off at a plateau, at ≈ 10% of its initial value after 400 fs. It then decays to zero on a much longer time scale, larger than 10 ps (see insert of Figure 4). This slow decay of the memory kernel compared to the VACF, which constitutes the main result of the present work, pinpoints the failure of the simple Langevin model, as the latter assumes an instantaneous relaxation of the solute-bath interaction compared to the dynamics of the solute (VACF). It may seem that the fast initial decay, which decreases the memory kernel by an order of magnitude, is sufficient to ensure a separation of time scales and the proper definition of a friction coefficient. This is however not the case, as the slower decay from this smaller value contributes significantly to the integral of K, which would define the friction if such a concept were applicable. Before analyzing in more detail the microscopic origin of this slow decay of the memory kernel, let us first compare the dynamics of confined Cs + with that in bulk water. Figure 5 compares the normalized VACF, FACF and memory kernel in the bulk and confined cases. The bulk results for the VACF and FACF are qualitatively similar to the confined case, even though quantitative differences are observed. For example, the decay of the VACF is slightly slower in the bulk and the negative minimum is only ≈ -0.2k B T /m. The initial decay of the FACF is almost identical to the confined case, but the subsequent oscillations are less pronounced. The impact of confinement is the largest on the memory kernel, which is the most important of the three correlation functions for the purpose of deriving continuous solvent models. The memory kernel, which follows the same initial sharp decay in the bulk and under extreme confinement, does not exhibit a plateau on this time scale and decays much more rapidly than in the confined case. Note that even in this case the contribution of the power-law decay (not visible on the figure) of the memory kernel to the integral which defines the friction coefficient is non negligible. In the bulk, this corresponds to the hydrodynamic backflow of solvent, which occurs on a time scale longer than the typical collision time with the solute with its neighbouring solvent molecules reflected in a negative minimum of the VACF for ≈ 0.5 ps.

Ion-water, ion-surface and ion-ion interactions

The sharper decay and oscillations of the VACF under extreme confinement suggests that the interaction of Cs + with water molecules is more repulsive in that case. The confined water molecules cannot transfer momentum to the rest of the fluid as easily as in the bulk and the solvent cage around the ion is stronger. In addition, the slower decay of the memory kernel also reflects the hindered momentum transfer to the rest of the fluid. In order to further quantify the role of ion-clay (IC), ion-water (IW) and ion-ion (II) interactions on the dynamics of the confined ion, we now analyze their correlation functions. Figure 6 reports the auto-and cross-correlation functions for these three specific forces, along and perpendicular to the clay surfaces.

In the direction along the surface, all these correlation functions do not show any significant variation over the picosecond time scale, with the notable exception of the ion-water autocorrelation function, which displays a minimum for ≈0. demonstrates that the evolution of the FACF and of the memory kernel discussed in the previous section is due to the ion-water interactions and more precisely to the oscillations of the Cs + ion in a solvent cage which otherwise does not evolve on this time scale. One can futher notice that the order of magnitude of all forces is similar: The initial value of the auto-correlation functions increases by a factor of 4 from ion-clay to ion-ion and ion-water interactions (IC-IC, II-II and IW-IW on Figure 6, respectively), indicating a two-fold increase in F 2 . In addition, the ion-clay and ion-water interactions are almost uncorrelated (small values of the cross-correlation function) , while the ion-ion interactions are anti-correlated with both the ion-water and ion-clay interactions (negative cross-correlation functions). Since the environment of the ion is frozen on this time scale, the ion-ion interaction and the parallel component of the ion-clay interaction can be considered constant and the anti-correlation only reflects the fact that the total force vanishes on average.

The force experienced by the confined ion in the direction perpendicular to the surfaces is, as expected, dominated by the ion-clay interaction, since the vertical confinement results from the steric repulsion by the clay surface. The perpendicular component of ion-ion interaction is negligible compared to the ion-clay and ion-water interactions, because the ion-ion vectors are approximately parallel to the surfaces. In that direction, only the ion-clay interaction evolves significantly over a few ps as the ion oscillates on this time scale, as indicated by the weakly damped oscillatory behaviour of the perpendicular component of this force. Note that it does not decay to zero on this time scale, reflecting the asymmetry of the instantaneous local environment of the ion near the middle of the interlayer. Thus the description of the system as purely twodimensional is only approximate. As for the parallel component, the anti-correlation between the ion-clay and ion-water interactions reflects the vanishing of the total force on average over this time scale.

Long-time dynamics

All these correlation functions decay to zero on a much longer time scale, as illustrated on Figure 7. Despite the decorrelation of the total force within 1 ps, it takes several ns for the ion-water, ion-clay and ion-ion interactions to decorrelate. The picture that emerges from this observation, together with the previous results, is that of a multi-scale process, with a short equilibration of the ion inside a local free energy minimum corresponding to a surface cavity and a solvent cage, and a much slower relaxation of this local environment and the corresponding motion of the ion to another local free energy minimum. Such slow relaxation is reminiscent of that observed in glasses or with site-to-site hopping in solids. The theoretical tools used to model the dynamics of such systems may thus prove useful in the description of the long-time dynamics of solutes under extreme confinement.

Conclusion

Using the Mori-Zwanzig formalism of projection operators, we have derived a Generalized Langevin Equation for the continuous solvent description of the dynamics of a confined solute. The solute evolves in the Potential of Mean Force under the effect of a random noise and a friction with a memory that both differ from that in a bulk solvent. As an example, we have considered the case of the Cs + ion confined in the interlayer of montmorillonite clays with a very low water content. We have computed the memory kernel from molecular dynamics simulations and have shown that it decays on a time scale longer that 10 ps, while the main features of the velocity autocorrelation function decay within a few ps. Thus there is no separation of time scales between the dynamics of the confined Cs + and that of the confined solvent, despite the much larger mass of the ion. This explains why the dynamics of this ion under extreme confinement cannot be described by a simple Langevin equation, as was previously attempted.

In the bulk, the slowest decay mode for the memory function, which is faster than in the confined case, corresponds to the hydrodynamic backflow of the solvent, while the shortest corresponds to the collision of the solute with its solvent cage. In order to understand the microscopic origin of the slow decay of the memory function under extreme confinement, we have further analyzed the relative contribution of the ion-surface, ion-solvent and ion-ion interactions to the dynamics. On the ps time scale, the only appreciable motion corresponds to the vibration of the ion in its neighbour cage, which involves the clay atoms in the direction perpendicular to the surface and water molecules in the direction parallel to the surfaces. This cage relaxes on a much longer (up to several ns) time scale, with the hopping of ions between local free energy minima, probably via the correlated motion of ions and solvent molecules. The resulting overall dynamics ressembles that of glasses or diffusion inside a solid by site-to-site hopping. This work opens a new perspective for the description of the dynamics of ions in water under extreme confinement, which could be applied in other materials of great scientific and practical relevance, such as zeolites [START_REF] Jeffroy | Understanding the Equilibrium Ion Exchange Properties in Faujasite Zeolite from Monte Carlo Simulations[END_REF] or carbon nanotubes [START_REF] Bonthuis | Theory and simulations of water flow through carbon nanotubes: prospects and pitfalls[END_REF]. Since in the bulk a Langevin description is relevant, it would be interesting to investigate wider pores, such as silica slit pores with a few water layers [START_REF] Argyris | Ion-Specific Effects under Confinement: The Role of Interfacial Water[END_REF], in order to determine where the cross-over between the extreme confinement and bulk-like regimes takes place.

Figure 1 .

 1 Figure 1. Cesium ions confined between montmorillonite clay surfaces. The interlayer porosity, seen from the side (left) and the top (right) also contains a monolayer of water. Cyan: Cs, yellow: Si, green: Al and Mg, red: O, white: H, grey: water molecules.

  , contains two layers of montmorillonite (in the figure, one of these layers is divided into two halves shown on the top and bottom) of 8 × 4 unit cells each of unit cell formula Cs 0.75 Si 8 (Al 3.25 Mg 0.75 )O 20 (OH) 4 . The 48 Cs + counterions, together with 288 water molecules, are located in the two interlayer porosities. The interlayer distance is 12.6 Å and the two layers are positioned so that the hexagonal cavities on their surface, visible in Figure 1, face each other. The lengths of the simulation box in the x and y directions parallel to the clay layers are 41.44 × 35.88 Å2 . For the aqueous Cs + ion, the cubic box of length 18.65 Å contains 1 ion and 215 water molecules. Periodic boundary conditions are applied in all directions.

Figure 3 .

 3 Figure 3. Two-dimensional Potential of Mean Force (PMF) in units of the thermal energy k B T for an interlayer cesium ion.

Figure 4 .

 4 Figure 4. Normalized autocorrelation functions for the velocity (VACF), force (FACF) and noise (memory function K), for an interlayer cesium ion. The insert shows the slower decay of the memory function compared to the VACF and FACF.

4. 3 Figure 5 .

 35 Figure 5. Normalized autocorrelation functions for the velocity (VACF), force (FACF) and noise (memory function K): comparison between a cesium ion in the interlayer of montomorillonite clay (lines without symbols) and in bulk water (lines with symbols).

Figure 6 .

 6 Figure 6. Specific force correlation functions (SFCF), for the three contributions to the force acting on the ion: ion-clay (IC), ion-water (IW) and ion-ion (II) interactions. Both auto-(IC-IC, IW-IW, II-II) and cross-(IC-IW, IW-II, II-IC) correlation functions are reported. Left Components parallel to the surface. Right Component perpendicular to the surface.

Figure 7 .

 7 Figure 7. Correlation functions of the different forces acting on the ion. The total force is separated into ion-clay (IC), ion-water (IW) and ion-ion (II) interactions. Both auto-(IC-IC, IW-IW, II-II) and cross-(IC-IW, IW-II, II-IC) correlation functions are reported, for the components parallel to the surfaces.

  1 ps followed by a plateau. This clearly

	2	1.5e+07					2	9e+06
	SFCF in (kJ/mol/Å)	0 -1.5e+07 -1e+07 -5e+06 0 5e+06 1e+07	0.2	0.4 t (ps) 0.6	0.8 IC-IC IW-IW II-II IC-IW IW-II II-IC	1	SFCF in (kJ/mol/Å)	0 -3e+06 0 3e+06 6e+06	1	2 t (ps) 3	4 IC-IC IW-IW II-II IC-IW IW-II II-IC	5

Acknowledgements

The authors thank Rodolphe Vuilleumier and Jean-François Dufrêche for useful discussions. AC acknowledges financial support from the Ecole Normale Supérieure (Paris). VM, MS, JPH, PT and BR acknowledge financial support from the Agence Nationale de la Recherche under grant ANR-09-SYSC-012.


Concerning the plausibility of the findings reported in 'Prenatal exposure to glycol ethers and cryptorchidism and hypospadias a nested case-control study' by Smet and Kelsey authors' response

Charline Warembourg, Luc Multigner, Ronan Garlantezec

Concerning the plausibility of the findings reported in 'Prenatal exposure to glycol ethers and cryptorchidism and hypospadias a nested case-control study' by Smet

Conflict of interest:

None to declare

We read with attention the letter from Smet and Kelsey discussing the interpretation of our findings. Smet and Kelsey 1 questioned the validity of our findings regarding the high detection rate of methoxyacetic acid (MAA) we observed in our population despite the restriction in use of its precursors in Europe. First, human biomonitoring studies on glycol ether (GE) exposure in the non-occupational context are very few but other recent studies have also reported a high detection rate of MAA in the European general population in 2007-2009: MAA was detected in 99% and 100% of urine samples with a median concentration of 0.05 mg/L and 0.11mg/L in a German 2 and a French population, [START_REF] Nisse | Caractérisation de l'exposition aux éthers de glycol d'un échantillon de population générale du Nord-Pas-de-Calais par biométrologie urinaire[END_REF] respectively. Second, regarding the sources of exposure to diethylene glycol monomethyl ether (DEGME) and triethylene glycol monomethyl ether (TEGME) as potential precursors of MAA detected in our study population, the period of recruitment (2002-2006) of the pregnant women included in our study should be considered.

The references given by Smet and Kelsey are not reflecting the usage at the time of the study. Indeed, DEGME is no longer used in consumer products since 2010 (at least without a concentration exceeding 0.1% by weight) but was previously used in paints and cleaning agents as assessed by the European chemicals agency [START_REF]Annex XVII to REACH -Conditions of restriction of DEGME[END_REF] and the French agency for environmental and occupational safety. [START_REF] Afsset | Les éthers de glycol -Synthèse des connaissances sur les expositions de la population générale en France[END_REF] Nowadays, as Smet et al said, TEGME is mainly use in hydraulic fluids but the reference they cited reported additional usages in consumer products such as in coatings, adhesives or inks, and that 'other release to the environment is likely to occur from […] indoor use (eg, machine wash liquids/detergents, automotive care products, paints and coating or adhesives, fragrances and air fresheners). [START_REF]Brief profile for TEGME[END_REF] Finally, we would like to highlight that each of these two GE (DEGME and TEGME) are manufactured or imported in 10000 to 100000 tons per annum in Europe. [START_REF]Brief profile for TEGME[END_REF][START_REF]Brief profile for DEGME[END_REF] Smet and Kelsey consider that the MAA we measured could not come from GE exposures; however, they did not provide any alternative exposure source that may explain the detection of MAA in urine. To our knowledge, there are no other compounds than GE as precursors of MAA.

As a second limitation, Smet and Kelsey indicate that the prevalence of male genital anomalies in our study is low. As we acknowledged in the discussion section of our article, the prevalence of male genital anomalies (especially cryptorchidism) is highly dependent of its definition (ie, the inclusion or exclusion of minor cases). [START_REF] Boisen | Difference in prevalence of congenital cryptorchidism in infants between two Nordic countries[END_REF] In our study, 1.9% of boys were diagnosed with cryptorchidism at birth (excluding cases with testis in high scrotal position) and 0.4% with hypospadias. Comparing with data published by Boisen et al (original reference used by Skakkebaek), the prevalence of cryptorchidism after exclusion of cases with testis in high scrotal position is similar to the one we reported: 1.2% in Finland and 2.1% in Denmark. For hypospadias, prevalence in the EDEN (Etude des déterminants pré et postnatals du développement et de la santé de l'enfant) and PELAGIE (Perturbateurs endocriniens: étude longitudinale sur les anomalies de la grossesse, l'infertilité et l'enfance) cohorts are in accordance with existing prevalence data. [START_REF] Springer | Worldwide prevalence of hypospadias[END_REF] In any case, it is known that the toxic effect of GE is caused by their alkox acid metabolites-in particular, MAA which is a known testicular toxicant-and ignoring the parent compound did not question our conclusion. Additional human biomonitoring studies would be important to assess current exposure level and evaluate the source of GE metabolites. Such studies will help to correctly perform a risk assessment evaluation and potentially provide recommendations to limit exposure of pregnant women. In addition, even if the European population is no longer being exposed to MAA, our results are of concern for non-European countries where precursors of MAA are not regulated.

and Kelsey: authors' response Reply Letter to Smet and Kelsey:

  Concerning the plausibility of the findings reported in "Prenatal exposure to glycol ethers and cryptorchidism and hypospadias: a nested case-control study", Warembourg et al, Occup Environ Med 2018;75:59-65 Warembourg C 1 , Multigner L 1 , and Garlantézec R 1,2 .

	1 Univ Rennes, Inserm, EHESP, Irset (Institut de recherche en santé, environnement et travail) -
	UMR_S 1085, F-35000, Rennes, France
	2 Service de Santé Publique et d'Epidemiologie, Centre Hospitalo-Universitaire de Rennes, F-
	35033, Rennes, France
	Corresponding author:
	Full name: Dr Charline Warembourg
	Postal address: UMR Inserm 1085 Irset -9 Avenue du Prof. Léon Bernard, 35000 Rennes,
	France
	Email: charline.warembourg@gmail.com


Benjamin Rotenberg

Olivier Bernard

Jean-Pierre Hansen

Salt-induced effective interactions and phase separation of an ultrasoft model of polyelectrolytes

We use a semi-grand canonical version of mean field density functional theory to determine the total effective interaction energy of a solution of penetrable polyions characterized by a gaussian charge distribution, in the presence of added salt. We then apply this effective representation of semi-flexible polyelectrolyte chains to investigate the possibility of a phase separation similar to that predicted earlier for charge-stabilized hard sphere colloids. Apart from the absence of a hard core repulsion, the effective pair potential is similar to the familiar DLVO potential between charged stabilized colloids, i.e. of the screened-Coulomb (Yukawa) form, but the effective valence of the polyions differs significantly from that of the DLVO pair potential, especially at high salt concentration. The existence of a well-defined closed-loop spinodal curve predicted by our mean-field calculation points to a phase separation betweeen solutions with high and low polyion concentrations under reasonable physical conditions. The salt concentration at the upper critical point is typically two orders of magnitude larger than in the case of hard core polyions, indicating that polyion penetrability appears to enhance the tendency towards phase separation.

Introduction

Linear polyelectrolytes are polymer chains carrying anionic or cationic groups which dissociate in water, releasing counterions of opposite charge. A typical example is polystyrene sulfonate which releases Na + or K + cations, leaving a negatively charged polyion. In the limit of low polyion concentrations, the counterions tend to move away from the polyion to maximize their entropy. The highly charged polyion then stretches into an elongated rod-like shape to minimize the electrostatic repulsion between the ionized segments (for reviews, see [START_REF] Barrat | Theory of Polyelectrolyte Solutions[END_REF][START_REF] Dobrynin | [END_REF]). At higher concentration, a fraction of the counterions will "condense" onto the charged polyelectrolyte (Manning condensation [3]), thus reducing the electrostatic repulsion between charged segmets, and inducing conformational changes by allowing the polyelectrolyte to fold into a quasi-spherical coil, or a "necklace" of such coils [START_REF] Barrat | Theory of Polyelectrolyte Solutions[END_REF][START_REF] Dobrynin | [END_REF]. Upon adding salt to the polyelectrolyte solution, the electrostatic interactions between segments are increasingly screened, thus enhancing the tendency towards coil-like conformations, as in the case of neutral polymers with short-range interactions between monomers.

Under such conditions it is not unreasonable to adopt a highly coarse-grained description of polyelectrolytes modelled as penetrable spherical objects characterized by a continuous gaussian charge distribution with a width of the order of the radius of gyration R g of the coil. Such a representation was recently introduced to investigate solutions of oppositely charged polyelectrolytes, within the so-called "ultrasoft primitive model" (UPM) [4][5][6][7]. In the present paper we adapt the model to solutions of equally charged polyelectrolytes in the presence of added salt. We use a classical density functional formulation, similar to that used in Ref. [8], to determine the equilibrium co-and counterion density distributions inside and around the penetrable polyelectrolyte coils, and to calculate the resulting effective pair interactions between "dressed" polyions.

The main objective of these calculations is two-fold. First of all we wish to compare the effective pair potential between penetrable polyions to the well-known DLVO pair potential between charged hard sphere colloids [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF]. In the latter case electric double-layers build up near the surface of the charged, inpenetrable colloid, while in the case of penetrable polyions considered here, the co-and counterions spread out over the interior as well as outside the polyelectrolyte coils, so that their local densities are no longer confined outside the polyions, and vary much more smoothly through space. Secondly, we will use the effective interaction between "dressed" polyions to calculate the free energy of the polyion solution, as a function of polyion and salt concentrations, within the random phase approximation (RPA), and explore the possibility of a phase separation between dilute and concentrated polyelectrolyte solutions [START_REF] Barrat | Theory of Polyelectrolyte Solutions[END_REF][START_REF] Dobrynin | [END_REF]. The theory described in the following sections of this paper closely follows a similar calculation for the case of dispersions of charged hard sphere colloids, which indeed predicts a re-entrant phase separation at low salt concentrations [START_REF] Van Roij | [END_REF][11][12]. The main interest of the present investigation is to analyse how polyion penetrability may affect such a phase separation.

A coarse-grained model of polyelectrolyte coils with added salt

We consider an aqueous solution of monodisperse polyelectrolyte chains of total charge Ze and monovalent point counterions of charge -e (assuming Z > 0 without loss of generality), in osmotic equilibrium with a salt reservoir of ideal monovalent co-and counterions. Within the usual "primitive model" of ionic solutions, the polar solvent reduces to a dielectric continuum of permittivity (electrostatic units will be used throughout). The polyelectrolyte solution is hence a three-component "mixture" of n 1 polyions (valence Z 1 = Z), n 2 counterions (valence Z 2 = -1) and n 3 coions (valence Z 3 = +1) per unit volume. Overall charge neutrality implies:

3 α=1 Z α n α = Zn 1 -n 2 + n 3 = 0 . (1)
In the coarse-grained representation [4,8], the polyelectrolyte coils are characterized by a "quenched" gaussian charge distribution of width σ 1 :

ρ1 (r) = Z 1 e 1 2πσ 2 1 3/2 e -r 2 2σ 2 1 , ( 2
)
where r is the distance from the centre of mass position R i of the ith coil (1 ≤ i ≤ N 1 ). A "quenched" (or "frozen") distribution is obviously an approximation since it implies that the internal structure (form factor) of the coil is unaffected by fluctuations of the local charge density of overlapping polyions or microions. The (

The electrostatic potentials generated by these distributions are:

ϕ α (r) = Z α e r erf r √ 2σ α ; 1 ≤ α ≤ 3 (4)
and the resulting bare pair potentials between ions are [8]:

v αβ (r) = ϕ α (r )ρ β (|r -r ) dr = Z α Z β e 2 r erf r √ 2σ αβ , (5)
where

σ 2 αβ = 1 2 (σ 2 α + σ 2 β
). Note that the electrosatic repulsion between two polyions and the attraction or repulsion between a polyion and a microion at full overlap v 11 (r = 0) = u 11 = Z 2 e 2 / √ π σ 11 and v 12 (r = 0) = -v 13 (r = 0) = -Ze 2 / √ π σ 12 are finite, thus preventing the "Coulomb collapse" of polyions and counterions.

The chemical potential µ 2 = µ 3 = µ s of the co-and counterions is determined by the salt concentration of the reservoir; if the reservoir ions are assumed to be ideal, their chemical potential µ s = k B T ln(n r s Λ 3 ), where n r s = n r 2 = n r 3 is the reservoir concentration of anion-cation pairs and Λ is the (irrelevant) thermal de Broglie wavelength of the ions.

Density functional theory of the effective interactions between polyions

The usual strategy for determining the effective interactions between polyions is to trace out the degrees of freedom of the microions for a given configuration {R i } of the former. Contrary to most previous work based on a canonical description of a closed system of polyions and microions, we adopt a semi-grand canonical point of view which incorporates the equilibrium with the salt reservoir [13,14]. The independent thermodynamic state variables are hence the total volume of the solution V , the temperature T , the total number N 1 of polyions (or equivalently n 1 = N 1 /V ) and the salt chemical potential µ s (or equivalently the reservoir salt concentration n r s ). Accordingly, the total effective interaction between polyions is [START_REF] Hansen | Theory of Simple Liquids[END_REF]:

V ef f 11 ({R i }) = V 11 ({R i }) + Ω({R i }) , (6)
where

V 11 ({R i }) = i j>i v 11 (|R i -R j |)
is the bare Coulomb interaction between all polyions (cf. Eq. 5) and Ω is the equilibrium grand potential of the inhomogeneous fluid of microions in the "external field" due to the polyions located at {R i }; henceforth the dependence of Ω on T and V will be understood. The grand potential is a functional of the local densities ρ 2 (r) and ρ 3 (r) of counterand coions:

Ω [ρ 2 (r), ρ 3 (r); µ s , {R i }] = F [ρ 2 (r), ρ 3 (r); {R i }] -µ s α=2,3 ρ α (r)dr , ( 7
)
where F is the free energy functional, which splits into ideal, mean-field and correlation contributions. The latter is neglected within the mean-field framework adopted in [8,[START_REF] Van Roij | [END_REF]11] as well as here. Moreover the ideal contribution to F is expanded to second order in the deviations ∆ρ α (r) = ρ α (r) -n α of the local densities of counterions (α = 2) and coions (α = 3) from their (yet unknown) mean values n 2 and n 3 . This approximation is adequate for slowly varying density profiles, a condition that is a priori better satisfied in the present case of penetrable polyions, compared to that of hard core polyions like charge-stabilized colloids examined in [START_REF] Van Roij | [END_REF][11][12]. The resulting approximate free energy functional reads [8,[START_REF] Van Roij | [END_REF]11]:

F [ρ 2 (r), ρ 3 (r); {R i }] = α=2,3 F id (n α ) + k B T ln(n α Λ 3 ) ∆ρ α (r)dr + k B T 2n α [∆ρ α (r)] 2 dr + e 2 ρ c (r) [ψ c (r) + 2ψ 1 (r)] dr , (8)
where

F id (n α ) = N α k B T ln(n α Λ 3 ) -1 , eρ c (r) = e [ρ 3 (r) -ρ 2 (r)] is the local mi- croion charge density, ψ c (r
) is the resulting electrostatic potential and ψ 1 (r) is the electrostatic potential due to the polyions at {R i }. ψ c satisfies Poisson's equation:

∇ 2 ψ c (r) = - 4πe ρ c (r) , (9)
while

ψ 1 (r) = N1 i=1 ϕ 1 (|r -R i |)
with ϕ 1 defined by Eq. ( 4). Substituting (8) in Eq. ( 7), minimization of Ω[ρ 2 , ρ 3 ] with respect to the two density profiles leads to the two coupled Euler-Lagrange equations:

k B T ln(n α Λ 3 ) + k B T n α ∆ρ α (r) + Ze [ψ c (r) + ψ 1 (r)] = µ s ; α = 2, 3 (10)
Two linear combinations of ρ 2 (r) and ρ 3 (r), namely ρ c (r) = ρ 3 (r) -ρ 2 (r) and ρ X (r) = x 2 ρ 3 (r) + x 3 ρ 2 (r), with x α = nα n (n = n 2 + n 3 ), decouple the two equations [START_REF] Van Roij | [END_REF] [START_REF] Van Roij | [END_REF]11], leading to:

∆ρ c (r) n = x 3 ln n r s n 3 -x 2 ln n r s n 2 -βe [ψ c (r) + ψ 1 (r)] (11)
and

ρ X (r) = n 2 n 3 n 2 + ln (n r s ) 2 n 2 n 3 , (12)
where β = 1/k B T and n r s is the reservoir salt concentration. ρ X (r) is uniform (independent of r). The normalization condition:

1 V ρ X (r)dr = 1 V [x 2 ρ 3 (r) + x 3 ρ 2 (r)] dr = 2n 2 n 3 n (13)
implies, by comparison with Eq. ( 12), the Donnan equilibrium relation:

n 2 n 3 = (n r s ) 2 = n r 2 n r 3 ( 14
)
so that ρ X (r) reduces to 2n2n3 n . Combining Eq. ( 12) with the electroneutrality relation (1) immediately leads to the following result for the counterion and coion concentrations in the polyelectrolyte solution:

n 2 = 1 2 n 2 1 Z 2 + 4(n r s ) 2 + n 1 Z ; n 3 = n 2 -n 1 Z . (15)
Combination of the Fourier transforms of Eqs. ( 9) and ( 11) and of ψ 1 (r) leads immediately to the following solution for the Fourier transform ρc (k) of the microion charge density ρ c (r):

ρc (k) = N1 i=1 ρci (k) = - N1 i=1 Zκ 2 D k 2 + κ 2 D e -k 2 σ 2 1 2 e ik•Ri , (16)
where

κ D = 4πne 2 kBT 1/2
is the Debye wavenumber associated with the microions. Inverse Fourier transformation leads to the microion charge density in real space:

ρ c (r) = N1 i=1 ρ ci (r) = - N1 i=1 Zκ 2 D 4π|r -R i | f 2 (|r -R i |) , (17)
where the function f 2 (s) is defined by:

f 2 (s) = 1 2 e κ 2 D σ 2 12 e -κDs erfc κ D σ 12 - s 2σ 12 -e κDs erfc κ D σ 12 + s 2σ 12 . (18)
ρ c (r) is thus the sum of "orbitals" centered on the N 1 polyions. The local electrostatic potential within the polyelectrolyte solution ψ(r) = ψ c (r) + ψ 1 (r) then follows directly via Poisson's equation:

ψ(r) = N1 i=1 Ze |r -R i | f 2 (|r -R i |) . (19)
Comparison of Eqs. ( 17) and (19) shows that:

ρ c (r) = -βenψ(r) . (20)
Note however that Poisson's equation determines the potential ψ(r) only within a constant, say ψ 0 . If ψ 0 is chosen to be zero in the reservoir, then comparison of Eqs. ( 20) and (11) shows that they are compatible provided a constant ψ 0 is added to ψ(r), namely:

ψ 0 = k B T e n 3 -n 2 n + ln n r s n 3 . ( 21
)
ψ 0 > 0 may be identifed with the Donnan potential; as expected ψ 0 → 0 as n 1 → 0.

The counterion and coion densities can finally be calculated from Eq. ( 17) and ρ X = 2n2n3 n via the inverse relations ρ 2 (r) = ρ X -x 2 ρ c (r) and ρ 3 (r) = ρ X +x 3 ρ c (r). Examples of the local densities ρ 2i (r) and ρ 3i (r) around an arbitrary polyion i, for two reservoir salt concentrations n r s , are shown in Figure 1. Their qualitative behaviour is as expected, with the counterion attracted inside the penetrable polyion, while coions are mostly expelled. At the higher salt concentration the profiles are nearly uniform. The equilibrium density profiles and corresponding potentials may now be inserted in the free energy functional (8) to calculate the equilibrium free energy of the microions as a function of polyion coordinates {R i }. The task is greatly simplified by the relation (20) which leads to partial cancellations between the ideal and mean field contributions to the free energy. Gathering results:

F [ρ 2 , ρ 3 ] = F id (n 2 ) + F id (n 3 ) - ZN 1 2 k B T ln n r s n 3 + e 2 ρ c (r)ψ 1 (r)dr . ( 22
)
The first three terms are independent of the polyion positions {R i } and are part of a structure-independent "volume term". The last, Coulomb term may be rewritten as:

e 2 ρ c (r)ψ 1 (r)dr = e i<j ρ ci (r)ψ 1j (r)dr + e 2 N 1 ρ ci (r)ψ 1i (r)dr . ( 23
)
The second term in Eq. ( 23) is structure-independent and contributes to the volume term. Gathering the results of Eqs. ( 24), ( 7), ( 17), ( 21), ( 22) and ( 23), we arrive at the following expression for the total effective interaction energy between the N 1 polyions:

V ef f 11 ({R i }) = i<j v ef f 11 (|R i -R j |) + U 0 . (24)
The effective pair potential is, setting r = |R i -R j |:

v ef f 11 (r) = v 11 (r) + e ρ ci (r )ψ 1j (|r -r|)dr = Z 2 e 2 r f 2 (r; σ 12 → σ 1 ) , (25)
where the function f 2 is defined in Eq. ( 18), with σ 12 replaced by σ 11 = σ 1 . In the February 12, 2014 Molecular Physics article limit of point polyions σ 1 → 0, f 2 (r; σ 12 → σ 1 ) simplifies to e -κDr , so that v ef f 11 (r) reduces to the familiar screened Coulomb (or Yukawa) form. For finite polyion size σ 1 , the asymptotic behaviour at long distance r is v ef f 11 (r) → Z 2 e 2 e κ 2 D σ 2 1 e -κ D r r [8]. Note that the "apparent valence" of the polyions is Ze κ 2 D σ 2 1 /2 , which differs from the DLVO apparent valence for hard sphere polyions, namely Ze κDσ1 /(1 + κ D σ 1 ), except for κ D σ 1 1. The full volume term U 0 in Eq. ( 24) reads:

U 0 = k B T N 2 ln n 2 n r s + N 3 + ZN 1 2 ln n 3 n r s -(N 2 + N 3 ) -N 1 Z 2 e 2 κ D 2 erfc (κ D σ 1 ) e κ 2 D σ 2 1 . (26)
Just as its bare counterpart v 11 (r), the effective pair potential ( 25) between polyions remains finite as r → 0:

v ef f 11 (r) r→0 ---→ Z 2 e 2 √ π σ 1 a(κ D σ 1 ) -b(κ D σ 1 ) r σ 1 2 + O(r 4 ) , ( 27
)
where a(x) = 1 -√ πxe x 2 erfc(x) and b(x) = (1 -2x 2 )/12 + √ πx 3 e x 2 erfc(x)/6. Note that v ef f 11 (r = 0) is strongly reduced in the high salt concentration regime κ D σ 1 > 1, where a(κ D σ 1 )

1, while at low salt concentration (weak screening; κ D σ 1 1) v ef f 11 (r) goes over to the bare pair potential at small r (a(κ D σ 1 ) → 1) [8].

Structure and thermodynamics of the fluid of "dressed polyions"

The microion degrees of freedom having been traced out, we are now left with a onecomponent system of N 1 "dressed" polyions whose total interaction energy is given by Eqns. ( 24)- (26). The pair structure and thermodynamics of this fluid of polyions can now be calculated using the usual tools of liquid state theory [START_REF] Hansen | Theory of Simple Liquids[END_REF]. In view of the absence of hard core in the smoothly vayring effective pair potential (25), the simplest approach is based on the random phase approximation (RPA), which has proved successful for systems of penetrable particles, like the gaussian core model [START_REF] Stillinger | [END_REF], provided the reduced amplitude βv(r = 0) is not too large [7,17]. In fact, the RPA becomes asymptotically exact in the high density limit of strongly overlapping particles [18,19].

The RPA assumes that the direct correlation function c(r) of the fluid of polyions is given by:

c(r) = -βv(r) = - Z 2 l B r f 2 (r; σ 1 ) (28)
where l B = βe 2 / is the familiar Bjerrum length and we have dropped the subscript "11" and the superscript "eff". The corresponding polyion-polyion structure factor follows from the Ornstein-Zernike relation [START_REF] Hansen | Theory of Simple Liquids[END_REF]:

S(k) = 1 + n 1 ĥ(k) = 1 1 -n 1 ĉ(k) = k 2 + κ 2 D k 2 + κ 2 D + κ 2 D1 e -k 2 σ 2 1 , (29)
where κ 2 D1 = 4πZ 2 l B n 1 . The excess free energy of the fluid of polyions due to the February 12, 2014 Molecular Physics article pair interactions may be calculated via Kirkwood's coupling constant integration method [20] for fixed values of n 1 , n r s and T :

βF ex 1 N 1 = n 1 2 1 0 dλ [1 + h(r; λ)] βv(r)4πr 2 dr = n 1 2 1 (2π) 3 1 0 dλ ∞ 0 ĥ(k; λ)βv(k)4πk 2 dk + n 1 2 βv(k = 0) (30)
Using again the RPA ansatz ĉ(k; λ) = -λβĉ(k) and the Orstein-Zernike relation to determine ĥ(k; λ) one arrives, after a trivial integration over the coupling constant λ, at the following expression:

βF ex 1 N 1 = κ 2 D1 2κ 2 D - κ 2 D1 8π 3/2 n 1 σ 1 + κ 2 D1 κ D 8πn 1 e κ 2 D σ 2 1 erfc(κ D σ 1 ) - 1 4π 2 n 1 ∞ 0 ln k 2 + κ 2 D k 2 + κ 2 D + κ 2 D1 e -k 2 σ 2 1 k 2 dk (31)
The total free energy F per polyion of the solution follows by adding the ideal term and the volume term U 0 defined in Eqn. (26). The last term in the latter exactly cancels the third term in (31). Defining the dimensionless free energy per unit volume ϕ = βF σ 3 1 /V and introducing dimensionless variables q = kσ 1 ; q D = κ D σ 1 ; q D1 = κ D1 σ 1 and n α ≡ n α σ 3 1 , the complete expression for ϕ reads:

ϕ = n 1 [ln n 1 -1] + n 1 q 2 D1 2q 2 D - q 2 D1 8π 3/2 - 1 4π 2 ∞ 0 ln q 2 + q 2 D q 2 + q 2 D + q 2 D1 e -q 2 q 2 dq + Zn 1 2 ln n r s n 3 -n 2 1 Z 2 + 4(n r s ) 2 (32)
We recall that n 2 and n 3 are determined by n 1 and n r s via Eqn. (14). For fixed values of T and V , ϕ is a function of n 1 and µ s (or n r s ); it is the thermodynamic potential associated with the semi-grand canonical description of the polyion-salt system adopted from the start. The chemical potential of the polyions follows from the thermodynamic relation:

βµ 1 = ∂ϕ(n 1 , µ s ) ∂n 1 µs . (33)
The total grand potential per unit volume of the solution, ω(µ 1 , µ s ), then follows from the Legendre transformation:

ω(µ 1 , µ s ) ≡ f (n 1 , µ s ) -n 1 µ 1 , (34)
where

f = F/V = k B T ϕ/σ 3 1 .
Remembering that the osmotic pressure P = -ω, we arrive at the following expression for the dimensionless osmotic pressure p = βP σ 3 1 : We now address the question of a possible demixing of the polyelectrolyte solution into a dilute ("gas") and a concentrated ("liquid") phase for fixed chemical potential µ s of the salt. Such a phase separation has been predicted to occur in dispersions of charged hard sphere colloids at very low salt concentrations [START_REF] Van Roij | [END_REF][11][12] and there is some experimental evidence for such a phase transition [21,22]. The important question is whether a similar phase separation occurs in the case of penetrable polyions, i.e. in the absence of a hard core repulsion. Phase separation occurs when a plot of the free energy function ϕ(n 1 , µ s ) versus polyion density n 1 at fixed salt chemical potential µ s exhibits a concave region, corresponding to metastable or unstable thermodynamic states; these are eliminated by the usual Maxwell double-tangent construction which determines the densities n g 1 and n l 1 of the coexisting "gas" and "liquid" phases. According to Eqn. (33), the slope of the double-tangent determines the common chemical potential µ 1 of the coexisting phases.

p = n 1 βµ 1 -ϕ(n 1 , µ s ) . (35)
In a first instance we have addressed the simpler problem of locating the spinodal line, i.e. the locus of points in the (n 1 , µ s ) plane where the curvature:

ξ(n 1 , µ s ) = ∂ 2 ϕ(n 1 , µ s ) ∂n 2 1 µs (36)
vanishes. The spinodal separates regions of positive and negative curvature of the free energy function, corresponding respectively to the stable or metastable states, and the thermodynamically unstable states. The spinodal and binodal (coexistence) curves meet at the critical point, or at an upper and a lower critical point in the case of closed spinodal curves corresponding to re-entrant phase behaviours (for example see references [START_REF] Van Roij | [END_REF]12,[23][24][25][26]). Fig. 2 shows the spinodal in the (n 1 , n r s ) plane, calculated for Z = 1000 and σ 1 = 70 nm. It is re-entrant, as in the case of charged hard-sphere colloids [START_REF] Van Roij | [END_REF]12], i.e. the spinodal forms a closed loop surrounding the region of thermodynamically 1 ≈ 0.00075 and a reservoir salt concentration n r s ≈ 2.2 mM. The latter value is about two orders of magnitude larger than typical values found in the hard core case [START_REF] Van Roij | [END_REF]11]. On the other hand the reduced polyion concentration n 1 σ 3 1 at the critical point is significantly lower than for hard core polyions, presumably as a consequence of polyion penetrability. The lower critical point occurs at extremely low values of n r s and n 1 (n r s ≈ 2 10 -11 M and n 1 ≈ 10 -11 ), which are physically irrelevant. Fig. 2 also shows the spinodal for σ 1 = 70 nm and the lower valence Z = 500 (corresponding to a weakly ionized polyelectrolyte). The critical reservoir salt concentration n r s is found to scale roughly like Z 2 , while the critical polyion density is practically unchanged.

Conclusion

We have used a semi-grand canonical version of mean field density functional theory to determine the total effective interaction energy of a solution of penetrable polyions characterized by a gaussian charge distribution, in the presence of added salt. We then applied this effective representation of semi-flexible polyelectrolyte chains to investigate the possibility of a phase separation similar to that predicted earlier for charge-stabilized hard sphere colloids. It is worth stressing that our ultrasoft polyelectrolyte model involves only purely Coulombic interactions between polyions and microions alike. Our main conclusions can be summarized as follows.

a) Apart from the absence of a hard core repulsion, the effective pair potential is similar to the familiar DLVO potential between charged stabilized colloids, i.e. of the screened-Coulomb (Yukawa) form, but the effective valence of the polyions differs significantly from that of the DLVO pair potential, especially at high salt concentration.

b) The existence of a well-defined closed-loop spinodal curve predicted by our mean-field calculation points to a phase separation betweeen solutions with high and low polyion concentrations under reasonable physical conditions. The salt concentration at the upper critical point is typically two orders of magnitude larger than in the case of hard core polyions, indicating that polyion penetrability appears to enhance the tendency towards phase separation.

A brief critique of the shortcomings of the model used and of the approximations made in the present work is in order. The model assumes a "quenched" charge distribution, which is clearly a gross over-simplification. In real polyelectrolyte solutions the internal structure and charge distribution of the polyions responds to variations of physical conditions, which induce significant changes of the polyelectrolyte conformations. Account must be taken of counterion "condensation" [3] and of the dependence of the electrostatic persistence length on the Debye screening length [START_REF] Barrat | Theory of Polyelectrolyte Solutions[END_REF][START_REF] Dobrynin | [END_REF]. Within the highly coarse-grained description adopted in this paper, at the very least the polyion charge Ze and size σ 1 should be made dependent on salt concentration.

As regards the approximations made in the present work, the RPA used in section 4 leads to the analytic expression (32) for the free energy which requires only a simple quadrature for numerical calculations. At very low salt concentrations the reduced value of the effective pair potential at full overlap βv 11 (r = 0) can be larger than 10 2 , which is too high to ensure the RPA to be quantitatively reliable. Although our predictions are expected to be at least qualitatively valid, since in the relevant range of low polyion concentrations considered here the probability of two polyions to overlap is very low, systematic corrections to the RPA will be included in future work which will focus on explicit calculations of the phase coexistence

February 12

 12 to be point ions (σ 2 = σ 3 = 0), such that: ρα (r) = Z α eδ(r) ; α = 2, 3 .

FebruaryFigure 1 .

 1 Figure1. Counterion (ρ 2 ) and coion (ρ 3 ) densities around a polyelectrolyte of radius σ 1 = 70 nm and charge Z = 1000e, at a volume fraction n 1 σ 3 1 = 0.1 in water at room temperature (βe 2 / σ 1 = 0.01) and two reservoir salt concentrations n r s = 10 -4 and 10 -2 mol.L -1 . In the former case, the deviation from unity of the long-distance limit reflects the Donnan effect.

February 12 5 .

 125 Polyion phase separation?

Figure 2 .

 2 Figure 2. Spinodal curve for a polyelectrolyte of radius σ 1 = 70 nm with charge Z = 500e (blue) or 1000e (red) in water at room temperature (βe 2 / σ 1 = 0.01).

February

  The upper critical point occurs at n 1 σ 3

February


domain: chim phys

text:
Alexandru Botan

Virginie Marry

Benjamin Rotenberg
email: benjamin.rotenberg@upmc.fr

Diffusion in bulk liquids: finite-size effects in anisotropic systems

We investigate systematically the effect of the cell size and shape on the diffusion properties in MD simulations. Specifically, we consider a bulk LJ fluid in orthorhombic cells with one length differing from the other two. The components of the diffusion tensor display complex variations as a function of the two independent lengths and may even become in some cases larger than the macroscopic limit for a cubic cell. These results can be perfectly explained by a purely hydrodynamic theory, which extends results obtained previously for the isotropic case. We provide the explicit expression of the diffusion tensor, including the effect of the finite size of the diffusing particle. The simulation results follow a simple scaling as a function of box size and aspect ratio and the corresponding scaling functions are determined numerically. These findings should have implications for the practically more relevant case of confined fluids.

Introduction

The diffusion of a solute or a solvent molecule arises from momentum transfer with the surrounding fluid, a process also at the origin of hydrodynamic flows around the particle -hence of hydrodynamic interactions between solutes. This fundamental correspondence is evident in the Stokes-Einstein relation, D ∞ = k B T /6πηa, relating the diffusion coefficient of a spherical particle of radius a to the solvent viscosity η (and the temperature T , with k B Boltzmann's constant), which can be derived from macroscopic hydrodynamics assuming stick boundary conditions at the surface.

One consequence of practical importance is that the computation of diffusion coefficients with Molecular Dynamics (MD) simulation suffers from finite-size effects arising from hydrodynamic interactions between periodic images [1,2]. For a cubic simulation box, the overall effect for a bulk fluid is a decrease in the diffusion coefficient. This decrease numerically follows down to almost molecular sizes the following scaling:

D(L) = D ∞ - ck B T 6πηL , (1)
with L the box length and c ≈ 2.837. The derivation of this result, already reported by Dünweg and Kremer [2] and by Yeh and Hummer [1] uses the Oseen tensor to compute the flow generated by the particle under a force F and Ewald summation to account for periodicity. They found that most of this effect is in fact due to the enforcement of a background force density -F/V , with V the volume of the system, ensuring overall momentum conservation. Higher order corrections accounting for the finite size of the solute can be obtained by considering the Rotne-Prager instead of the Oseen tensor [1,3,4]. Most reported simulation data, even since the publication of these results, neglect this contribution and therefore underestimate the diffusion coefficient. Nevertheless, this scaling law has been verified numerically for a number of systems, including several popular water models [5][6][7], ionic liquids [8] or more complex fluids such as solutions of star polymers [9] and even in higher dimensional spaces [10]. This scaling then allows extrapolating the solvent diffusion coefficient to the macroscopic L → ∞ limit, or the correction to be applied for solutes [11,12]. In addition, the slope provides an estimate of the viscosity -a strategy that was used in particular to determine the viscosity of water in ab initio MD simulations [13].

More recently, it was observed that for anistropic simulation boxes, the diffusion tensor becomes anisotropic even for bulk fluids [6]. Such finite-size artefacts may become problematic for the study of confined fluids, a very active field of research due to the numerous applications in a wide variety of contexts (e.g. micro-and nanofluidics, nanoporous materials, etc), for which anisotropic simulation boxes are the rule rather than the exception. Note that in the case of a slit pore, one dimension is indeed finite and only the periodicity in the other directions results in spurious size effects. For an anisotropic bulk water system, Rozmanov and Kusalik found that replacing L by V -1/3 in Eq. 1 is not sufficient to describe the components of the diffusion tensor. Rather, they proposed for the considered geometry to use instead the square root of the cross-sectional area corresponding to each direction [6].

Here we investigate systematically the effect of simulation box anisotropy on the diffusion tensor. We show that hydrodynamic effects can explain the simulation results, including the fact that some components may exceed the extrapolated value for a cubic box in the infinite size limit. To that end, we follow the strategy of Ref. [1] and consider explicitely (a) the Rotne-Prager tensor instead of Oseen, to account for the finite size of the diffusion particle (b) all components of the diffusion tensor instead of the isotropic part. This is presented in Section 2, together with the simulation methods. Section 3 briefly reports the results for isotropic boxes. The main results of the present work, dealing with anisotropic boxes, are finally presented in Section 4.

Hydrodynamic theory and simulation methods

Hydrodynamics

The effect of hydrodynamics on the diffusion tensor in periodic systems, as well as explicit expressions for its isotropic component and for cubic boxes has been given by Yeh and Hummer [1]. For the sake of completeness, we first recall the main arguments of their derivation and then extend it to provide the expressions for the full tensor, including the effect of the finite size of the particle via the Rotne-Prager tensor. The diffusion tensor under periodic boundary conditions is given by:

D a PBC = D ∞ 1 + k B T lim r→0 [T PBC (r) -T ∞ (r)] , (2)
where T PBC and T ∞ denote the mobility tensors under periodic and unbounded conditions, respectively, and 1 is the identity matrix. These tensors relate the velocity of the fluid at a position r around the particle when the latter is submitted to an applied force F, while a compensating background force -F/V is distributed uniformly over the fluid, as v(r) = T(r) • F. For point particles, the mobility tensor in the unbounded case is the Oseen tensor:

T 0 ∞ (r) = 1 8πηr 1 + rr r 2 , (3)
where the superscript refers to the vanishing radius a = 0, while for the periodic case it is more conveniently written in reciprocal space as:

T 0 PBC (r) = k =0 e -ik•r ηk 2 V 1 - kk k 2 , (4)
where the sum extends over the reciprocal lattice vectors k, excluding k = 0. For spheres of radius a = 0, the Oseen tensor is replaced by the Rotne-Prager tensor [3] and the extension of Eq. 4 has been given by Hasimoto [4] and Beenakker [14]. The diffusion tensor can then be written as an Ewald sum by introducing a real number ξ > 0:

D PBC = D ∞ 1 + k B T ξ 6πη -6 + 40a 2 ξ 2 3 1 √ π 1 + M r (a) + 6π V ξ 3 M k (a) (5)
where we have introduced two tensors:

M r (a) = s =0 f 1 (s)1 + f 2 (s) ss s 2 (6)
f 1 (s) = 3 4 + a 2 ξ 2 2s 2 erfc s s + - 9 2 + 14a 2 ξ 2 + 3s 2 -20a 2 ξ 2 s 2 + 4a 2 ξ 2 s 4 + a 2 ξ 2 s 2 e -s 2 √ π f 2 (s) = 3 4 - 3a 2 ξ 2 2s 2 erfc s s + 3 2 -2a 2 ξ 2 -3s 2 + 16a 2 ξ 2 s 2 -4a 2 ξ 2 s 4 - 3a 2 ξ 2 s 2 e -s 2 √ π and M k (a) = t =0 e -t 2 /4 t 2 1 + t 2 4 + t 4 8 1 - a 2 ξ 2 t 2 3 1 - tt t 2 . (7)
The two sums extend over real and reciprocal space lattice vectors, respectively. For an orthorhombic box with lengths L x , L y and L z , to which we will restrict ourselves in the following, this reads:

s = n x ξL x e x +
with n x,y,z and p x,y,z integers and where s = ||s|| and t = ||t|| denote the norm of these vectors. Eqs. 5-9 provide the extension of Eq. 1 to finite size solutes and to anisotropic simulation boxes, where the components of the diffusion tensor may differ. Eq. 1 is recovered in the a → 0 and taking the isotropic part ( 13 Tr D PBC ). In the more general case we were not able to compute the Ewald sums analytically, however they can be easily determined numerically provided a suitable choice of ξ and number of terms in the sums 6 and 7. In the following, results are obtained by choosing ξ = 6.4/ min α L α , n x,y,z ∈ [-n max , n max ] with n max = 1 and p x,y,z ∈ [-p max , p max ] with p max = 2ξ min α L α + 1, which are sufficient to ensure good convergence. For a systematic study of the effect of anisotropy, we consider for simplicity orthorhombic systems where two lengths are identical. Without loss of generality, we note L x = L y = L and L z = L ⊥ (see Figure 1). In that case, the diffusion tensor is diagonal with D xx = D yy = D and D zz = D ⊥ . We simulate a Lennard-Jones fluid under the same physical conditions as the one considered in Ref. [1], namely a reduced density ρ * = ρσ 3 = 0. Table 1. Size of the simulated systems (see Figure 1), in units of the LJ diameter σ. The first line corresponds to cubic cells.

Simulation details

Molecular simulations are performed in the microcanonical ensemble for 10 4 t * , with t * = σ m/ the Lennard-Jones time unit, using a time step of 10 -3 t * , after an equlibration period of 10 2 t * in the canonical ensemble. For the latter, a Nose-Hoover thermostat with a time constant of 0.1t * is used. A cut-off of 2.5σ is used to compute the Lennard-Jones interactions. For each system, 8 independent runs (16 for the smallest system L = L ⊥ = 5) with different initial conditions are performed. For each simulation, the mean-square displacement in each direction is computed and the components of the diffusion tensor are determined from the slope for t < 4.10 3 t * . Reported results are averages over these independent runs, with errorbars computed using the block averaging method [15].

Since the hydrodynamic results involve the viscosity of the fluid, we also compute this quantity from equilibrium simulations using the Green-Kubo expression: (10) involving the auto-correlation function (ACF) of the off-diagonal components of the stress tensor σ αβ . The reported results for the viscosity correspond to the L = L ⊥ = 15 system. As previously noted in the literature [1], we not observe significant size effect on the viscosity. All simulations are performed using the LAMMPS simulation package [16].

η = V k B T ∞ 0 σ αβ (t)σ αβ (0) dt

Isotropic boxes

Results for isotropic boxes (L = L ⊥ = L) are reported in Figure 2. As expected, they follow the linear dependence in 1/L of Eq. 1). The value of the viscosity, η = 1.28 ± 0.03, is in good agreement with Ref. [1] and describes very well the slope of D vs 1/L. The extrapolated value for infinite box size, D ∞ = 0.312 ± 0.005 is also consistent with Ref. [1]. These values of η and D ∞ are then used in the following for the hydrodynamic calculations in the anistropic cases. For box sizes comparable to the LJ diameter, one observes deviations from Eq. 1, as expected for this limit of point particles. Nevertheless the effect of the finite particle size is moderate, as can also been seen in Figure 2 which reports the exact numerical result for a size a = σ (the relative correction is of order a 2 /L 2 , as shown in Ref. [1]). For L ≥ 10σ, this effect can be neglected, as we will do in the following. We now turn to the more general case of anisotropic boxes and consider first the evolution of the components of the diffusion tensor D and D ⊥ for fixed values of L , as a function of the perpendicular dimension L ⊥ . Figure 3 reports results for L = 10 and 40σ. One first observes that the x and y components are equal, and that D and D ⊥ differ, as expected from the symmetry of the simulation cell. In addition, none of these components follow the simple scaling (Eq. 1) of the isotropic case. D and D ⊥ may decrease or increase, simultaneously or in opposite directions, depending on the size and aspect ratio. Both curves intersect as expected when L ⊥ = L . A striking observation is that both components may exceed the value D ∞ for an infinite cubic cell, whereas hydrodynamic effects can only reduce the diffusion coefficient for isotropic case. Finally, we note that the hydrodynamic calculations for point particles (a = 0 in Eqs. 5-7) perfectly reproduce the simulation results. This indicates that although the anisotropic case is more complex than the isotropic one, all the finite size effects are due to hydrodynamics, at least in this regime where L and L ⊥ are larger than 5σ. All the above discussion also applies to the evolution of D and D ⊥ for fixed values of L ⊥ , as a function of the perpendicular dimension L , as can be seen for L = 10 and 40 in Figure 3.

Scaling with size and aspect ratio

It may seem surprising that some components exceed in some cases the value for an infinite cubic cell and even seem to diverge e.g. as L ⊥ → ∞ for L = 10σ (Figure 3). One might expect, for example, to recover D ∞ for an infinitely large system. For example, comparing the L ⊥ → ∞ (resp. L → ∞) for L = 40σ (resp. L ⊥ ) to the same limit for 10σ, one already sees that both D ⊥ and D approach D ∞ . However, for an anisotropic system the proper macroscopic limit should be taken at fixed aspect ratio α = L ⊥ /L . Examination of Eqs. 5-7 shows that both component of the diffusion tensor can be rewritten as:

D ,⊥ = D ∞ + k B T 6πηL g ,⊥ ξL , L ⊥ L , a L , (11)
where the two functions g ,⊥ can be deduced straitforwardly from the above expressions. Noting that the Ewald summation technique is independent of the choice of ξ (when both sums extend over the complete real and reciprocal space), these two functions are in fact independent of ξL . In addition, we have seen in the previous section that molecular simulation results are perfectly reproduced, for the cases considered here (a/L ∼ σ/L < 1/5), by the hydrodynamic result for point particles (a = 0). In this limit, we thus obtain the following scaling:

D ,⊥ = D ∞ + k B T 6πηL h ,⊥ L ⊥ L , (12)
where the two functions h ,⊥ now depend only on the aspect ratio α = L ⊥ /L . Both functions can be determined numerically and the isotropic case already provides h ,⊥ (1) = -c ≈ -2.837. Figure 5 reports the scaling functions h ,⊥ = (D ,⊥ -D ∞ )/(k B T /6πηL ) as a function of the aspect ratio for all simulated systems. They all collapse on a master curve (one for each component), which is perfectly described by the hydrodynamic result of Eqs. 5-7.

This figure not only confirms the hydrodynamic origin of the non-trivial finitesize effect on the diffusion tensor in anisotropic simulation cells, but also provides the evolution of the scaling functions h ,⊥ . For the perpendicular component, h ⊥ (red curve in Figure 5) decreases monotonously and plateaus for large aspect ratios (thin rod), with lim α→∞ h ⊥ ≈ -2.925, compared to c ≈ -2.873 for α = 1 (cubic). For small aspect ratios (flat cell), h ⊥ diverges approximately as α -2 . The finitesize correction is positive (D ⊥ > D ∞ ) for α < α 0 ⊥ ≈ 0.396, while it is negative otherwise. As for the parallel component, h (blue curve in Figure 5) decreases for 1) are reported (symbols) together with the hydrodynamic result Eqs. 5-7. Note the logarithmic scale on the x-axis.

α < α m ≈ 0.396 down to h (α m ) ≈ -3.788 and increases for α > α m . It diverges approximately as α -2 for α → 0 and as α for α → ∞. Finally, the finite-size correction is positive (D > D ∞ ) for α < α 1 ≈ 0.146 or α > α 2 ≈ 2.794 and negative otherwise. The fact that the minimum of h is obtained for a value α m approximately equal to the zero α 0 ⊥ of h ⊥ is striking, but we were not able to demonstrate the equality between these two particular points.

Conclusion

We have investigated systematically the effect of the cell size and shape on the diffusion properties in MD simulations. Specifically, we considered a bulk LJ fluid in orthorhombic cells with one length differing from the other two. We showed that the components of the diffusion tensor display complex variations as a function of the two independent lengths and may even become in some cases larger than the macroscopic limit for a cubic cell (for which the finite size effect is always a decrease in the diffusion coefficient). These variations are perfectly described by a purely hydrodynamic theory, which extends results obtained previously for the anisotropic case. In particular, we provided the explicit expression of the diffusion tensor including the effect of the finite size of the diffusing particle. The latter turns out to be small in the present case, but may be more important for solutes larger than solvent molecules in small simulation cells. Using this hydrodynamic description, we have further shown that the simulation results follow a simple scaling as a function of box size and aspect ratio and the corresponding scaling functions have been determined numerically. These findings should have implications for confined fluids (see e.g. Ref. [17]): In that case hydrodynamic flows in the directions parallel to the confining walls may result in spurious size effects on the computed diffusion coefficients. The approach followed in the present work should allow to adress this more practically important case.

Figure 1 .

 1 Figure 1. The simulated systems consist of a LJ fluid in orthorhombic cells with one length different from the other two. Both elongated (L ⊥ > L , as shown) and flat (L ⊥ < L ) are considered.

  7 and a reduced temperature T * = k B T / = 2.75, with σ and the Lennard-Jones diameter and energy, in the supercritical fluid region. The simulated systems, illustrated in Figure1, cover a range of sizes and aspect ratios from flat to elongated. Specifically, the box lengths corresponding to the results reported below are summarized in table1, with particle numbers ranging from 87 to 87500.

Figure 2 .

 2 Figure 2. Diffusion coefficient for isotropic systems (L = L ⊥ = L). Simulation results (blue symbols) are compared to the hydrodynamic calculations for a point particle (Eq. 1, black solid line) and with a size equal to the LJ diameter (a = σ in Eqs. 5-7, red dashed line). The inset shows the integral defining the viscosity in the Green-Kubo relation Eq. 10 (with errorbars in light grey). The corresponding value of the viscosity is used for the hydrodynamic calculations and the only fitted parameter is the extrapolation to infinite box size D∞ (indicated by a green square), used in the following. All quantities are given in LJ units.

Figure 3 .

 3 Figure 3. Diagonal components of the diffusion tensor, for Lx = Ly = L = 10 and 40, as a function of 1/Lz = 1/L ⊥ . Simulation results (symbols) are compared to the hydrodynamic result for a point particle (a = 0 in Eqs. 5-7, lines). The diffusion coefficient extrapolated for an infinite isotropic box size, D∞, is also indicated (green square). All quantities are in LJ units.

Figure 4 .

 4 Figure 4. Diagonal components of the diffusion tensor, for Lz = L ⊥ = 10 and 40, as a function of 1/Lx = 1/Ly = 1/L . Simulation results (symbols) are compared to the hydrodynamic result for a point particle (a = 0 in Eqs. 5-7, lines). The diffusion coefficient extrapolated for an infinite isotropic box size, D∞, is also indicated (green square). All quantities are in LJ units.

Figure 5 .

 5 Figure 5. Scaling functions h ,⊥ = (D ,⊥ -D∞)/(k B T /6πηL ) as a function of the aspect ratio α = L ⊥ /L . Simulation results for all systems (see Table1) are reported (symbols) together with the hydrodynamic result Eqs. 5-7. Note the logarithmic scale on the x-axis.

  n y ξL y e y + n z ξL z e z

							(8)
	t = p x	2π ξL x	e x + p y	2π ξL y	e y + p z	2π ξL z	e z

Acknowledgements

The authors would like express their gratitude to Jean-Pierre Hansen for his guidance and their pleasure to enjoy him as a colleague at UPMC. BR is particularly indebted to Jean-Pierre for his continuing scientific inspiration and friendship, from his first interactions with him as a visiting student in Cambridge to a shared office in Paris.