Collective density fluctuations and dynamics of ions in water studied by the interaction-site model of liquids

The collective excitations in water are studied based on the interactionsite model of liquids. Three collective modes, extracted from a generalized Langevin equation combined with the RISM theory of liquids, are identified as an acoustic mode and two optical modes. The drag force exerted on ions in water is described in terms of the response of these solvent collective excitations to the perturbation of ions. The ion-size dependence of the drag force, which has been a central issue in physical chemistry for long time, is studied in molecular detail based on the novel approach.


Introduction
Study of dynamics in liquid water has scored considerable progress in the past few decades.Especially, the molecular simulation method has made great contribution revealing unique nature of the dynamics, which has its origin in the structural feature characterized by the hydrogen-bond network.
Recently, two interesting papers have appeared, which concern collective excitations in water.One of those is related to the acoustic mode.In general, if the characteristic frequency ω of density fluctuations in liquid is plotted against the wavenumber k, a linear relation is obtained in the small wavenumber limit, namely lim k→0 ω = ck, where c is the sound velocity.In the case of water, the sound velocity is approximately 1500 m/sec.However, it has been disclosed by Teixeira and coworkers based on the neutron scattering measurement that there is a 'new' sound mode which propagates with velocity twice as fast as the 'ordinary' sound [1].The sound velocity is very close to that of ice.The existence of a similar sound mode with a large velocity has been predicted by Rahman and Stillinger from the molecular dynamics simulation [2].These authors related the 'new' sound mode to the dynamics of the hydrogen-bond network.On the other hand, such enhance in sound velocity with increasing wavenumber is a phenomenon normally observed in liquids, even in argon, which is known as the 'positive visco-elastic effect ' [3].Therefore, it is hard to draw a definite conclusion at this moment with respect to which picture is closer to the reality.An interesting question to be asked in this regard is how the crossover occurs between the two sound modes when the wavenumber is varied.More recently, an Italian group published an X-ray scattering result, showing that the crossover takes place in continuous fashion [4].It is not so easy for molecular dynamics simulation to study such a problem, since there is a serious limitation concerning the system size, or number of molecules to be simulated, which determines the lower bound of the wavenumber.It would be more appropriate to rely on the non-equilibrium statistical mechanics in this respect.
The other topic which has stimulated the field is collective excitations corresponding to the "optical mode."All other modes in which the characteristic frequency ω does not vanish in the k → 0 limit are collectively referred to as optical modes.Vallauri and coworkers extracted dispersion relations corresponding to the optical modes from the analysis of the dynamic structure factor S(k, ω) of water [5].The optical modes correspond apparently to the rotational motion of water, since the model for water molecules is a rigid body and the vibrational excitation is entirely prohibited.

Theory of liquid dynamics and collective excitations in water
The present study is based on the generalized Langevin equation which has played a central role in the non-equilibrium statistical mechanics."Generalization" of the Langevin equation, often mentioned in literature, features a phenomenological replacement of the friction and random forces by their counter parts accounting for the non-Markovian nature of dynamics: for example, "colored noise" instead of "white noise" for the random force.Such a description, however, does not meet the requirements of chemistry, since chemical nature of environmental solvent around a Brownian particle is too complicated to be explained by such phenomenological parameters involved in the theory as the relaxation time of memory.There has been a considerable effort devoted to the "generalization" of the Langevin theory based on the projection-operator method, which relates the friction and random forces with microscopic structure of liquids in terms of the density pair-correlation functions [6].Such methods enable us to treat liquid dynamics at a molecular level at least for simple liquids consisting of spherical particles.However, the problem is not so simple for molecular liquids which involve not only the translational motion but also rotational motion of molecules.A straightforward solution to the problem is to represent molecular center and orientation by six coordinates, and such a description has been developed traditionally in literatures.The simplest example is the Debye-Fick diffusion equation, which has been employed in the NMR relaxation theory.The method, however, has obvious difficulty when it is applied to chemical processes.Firstly, using angular coordinates for rotational motions makes mathematical formula too complicated even for small molecules.Secondly, such a treatment requires an explicit description for coupling of rotational and translational degrees of freedom, which is another nontrivial problem.Furthermore, it will be impossible to apply the method to a chemical reaction in which molecules in concern change their composition and/or geometry along a reaction path.The consideration stated above made us construct a new model for dynamics of molecular liquids.What if we regard dynamics of molecules as "translational" motions of constituent "atoms", which are correlated with each other?The Site-Site Smoluchowski-Vlasov (SSSV) theory has been proposed based on such a consideration, in the most primitive level of overdamped regime [7].The model is a natural extension of the interaction-site model in the equilibrium statistical mechanics to non-equilibrium process, and the theory can be realized by combining the RISM theory with the generalized Langevin theory.The model featuring a correlated translational motion of atoms has its own disadvantage: the rotational motion of molecules is not explicitly described.The problem has been solved by Chong and Hirata by extracting collective modes of the density fluctuation from the site-site density correlation functions.In their recent study for dynamics of a simple dipolar liquid based on the interaction-site model [8], they have succeeded to abstract the collective excitations in liquids, which can be identified as optical and acoustic modes, by diagonalizing the collective frequency matrix appearing in the generalized Langevin equation.The two modes are related essentially to the rotational and translational motions of molecules.In what follows, we describe briefly the method applied to water.
Let us begin with the definition of the site-density fluctuations, where r α i (t) denotes the position of α atom in ith molecule at time t, and N is the total number of molecules.Here, α refers to oxygen atom (denoted as O) or one of two hydrogen atoms (H 1 and H 2 ) constituting a water molecule.The site-site density correlation function is defined by the initial value of which is the site-site static structure factor, where K(k, t) represents the memory function, and the matrix notation is used for site-site correlation functions such as F αβ (k, t).ω n k denotes the normalized nth frequency-moment matrix of S(k, ω) defined by the relations, If the damping term is absent, equation ( 3) is formally identical to that of the harmonic oscillator, q + ω 2 q = 0.The eigen modes of the collective density fluctuations can be obtained by diagonalizing the matrix ω 2 k , whose explicit expression for a three-site water model can be found elsewhere [9].The eigen frequencies obtained by diagonalizing ω 2 k for water are plotted against the wavenumber k in figure 1a.There are three eigen frequencies, one of which vanishes in the k → 0 limit.By definition, the mode can be identified as an 'acoustic' mode.The asymptotic form of the eigen frequency of the mode becomes in the non-damping case [9] as follows: where M is the total mass of a molecule, and we have noticed that all the sitesite static structure factors coincide in the k → 0 limit, denoted as χ(k = 0).The ω 2 acou (k → 0) coincides with that of the ordinary sound mode propagating with the isothermal sound velocity (k B T /Mχ(k = 0)) 1/2 and the wavenumber k [10].The fact that the above expression includes the total mass of a molecule suggests the mode is related to the translational motion of molecules.
The other two eigen frequencies, on the other hand, do not vanish in the k → 0 limit.By an analogy to the solid state physics, the other modes can be identified as 'optical' ones.These modes are concerned with the relative motions of atoms, whose characteristic frequencies do not vanish even in the k → 0 limit.The limiting behaviours of the two optical modes, to be called OM-I (the optical mode I) and OM-II (the optical mode II), are respectively given by [9] where I x , I y and I z are the moments of inertia around the principal x-, y-and z-axes, and z O and z H are the z-coordinates of O and H atoms of a molecule in the body-fixed coordinate (see figure 1 of [9]).The optical modes are related to the rotational motion of molecules as can be inferred from the appearance of the moments of inertia in the above expressions.OM-I has a the second moment of the density correlation function, χ ′′ (k = 0), in its expression, which is closely related to the dielectric constant of liquid [11,12].It also indicates that OM-I has collective character.On the other hand, OM-II is essentially a single-molecule mode since the expression is not associated with any collective density correlation function.
The contributions from each atom to the mode can be extracted in the following way.Diagonalizing the matrix ω 2 k corresponds to turning the description of the system in terms of the density fluctuations ρ α (k) to the one in terms of their linear combination: where x O (k), x H 1 (k) and x H 2 (k) are the components of the eigen vector corresponding to the mode.Thus, by analyzing the sign and magnitude of x α (k)'s, it is possible to obtain the information on how each atom contributes to the mode.Plotted in figure 1b to 1d are the contributions of each atom to the acoustic and two optical modes, respectively.It is seen from figure 1b for the acous- holds well in the small-k region, which is consistent with the fact that the sound mode stems from the center-of-mass translational motion of the molecules, i.e., each atom in the molecule equally contributes to this mode.On the other hand, as can be seen from figure 1c, OM-I is governed by the lighter hydrogen atoms over the entire wavelength range, because the rotational motion of a molecule is dominated by the motion of the lighter atoms which are located further from the center of mass. Figure 1d shows that OM-II is just related to the rotational motion in which two hydrogen atoms move in the out-ofphase fashion with oxygen atom fixed.More detailed analysis of the eigen vectors revealed that OM-I and OM-II are associated with the pitch and roll librational motions of the molecules, respectively (see figure 4 of [9]).Since OM-II turns out not to contribute to the ion dynamics to be described in the next section, it shall not be considered any more in the present article.
These resonances determined by ω 2 k are shifted and damped by the memory functions K αβ (k, t).An approximation scheme has been developed for K αβ (k, t) which enables one to calculate the site-site density correlation functions F αβ (k, t) and their spectra [8].This comprises the assumption of the exponential form for the memory functions and the extension of Lovesey's method for determining relaxation times appearing in the assumed form for the memory functions.Figures 2 and 3 show the results based on this theory along with the molecular-dynamics (MD) simulation data.The results are reported in the form of longitudinal-current spectra, which give the spectra of the collective excitations in liquids, defined by Furthermore, linear combinations of the form rather than C L,αβ (k, ω) themselves are plotted in the figures since one can separately discuss the acoustic and optical modes based on these combinations [13].Here, C L,MM (k, ω) and C L,ZZ (k, ω) denote the longitudinal-current spectra of the mass and charge, respectively, which can be obtained by setting c α = m α (mass of atom α) or c α = q α (partial charge of atom α) in equation ( 11).It is seen from figure 2 that the theoretical results for the acoustic-mode spectra are in fair agreement with MD simulation results; especially, the peak positions of C L,MM (k, ω) are well reproduced.(The discernible lower frequency peak in MD simulation result for n = 4 in figure 2 can be ascribed to a single-molecule excitation, and not to the collective acoustic one, as exhibited by Miura [14].)Concerning the width of the spectra, the agreement is not so satisfactory, and the use of the recently developed mode-coupling theory for molecular liquids [15,16] may improve the results in this regard.Compared in figure 3 are the theoretical and MD simulation results for the longitudinal current spectra of the optical dynamics, C L,ZZ (k, ω).As can be seen from figure 3, the theoretical results for the peak frequencies are considerably lower than those of MD simulation.However, the overall shape of the spectra is well reproduced by our theory as can be checked by shifting the theoretical results horizontally such that their peak positions coincide with those of the simulation results.It was discussed that the disagreement in the peak positions of C L,ZZ (k, ω) is largely due to the disagreement between the theoretical and simulation results for the fourth frequency moment ω 4 k [13].Thus, although some quantitative disagreements can be seen between the theoretical and MD simulation results, the qualitative nature of the collective density excitations in liquid water is well captured by the present theory.Some other features on collective excitations in water, such as the high-frequency sound velocity and the dispersion relations, have also been discussed in [9,13], showing that all the essential features, reported previously by neutron-scattering experiments [1], MD simulations [5,14,[17][18][19] and dielectric theories [20,21] are well reproduced.

Dynamics of solvated ions in water
Dynamics of ions in polar liquids has attracted scientists' attention for a long time in many different fields including physics, chemistry and biology.Main interests have been focused upon the puzzling behaviour exhibited in the ion-size dependence of the friction coefficient, which decreases with increasing ion-radii as opposed to the prediction of the hydrodynamic Stokes law.Two classical models which are based on entirely different physics have been practiced to explain the paradoxical behaviour.The first of those is the so-called 'solvent-berg' model in which an ion makes a 'complex' with solvent molecules, and undergoes Brownian motion with the 'effective Stokes radius' of the molecular complex.The effective radius decreases with increasing radius of a bare ion, so does the Stokes friction, because the electrostatic attraction between ion and solvent decreases with increasing ionic radius.The dependence of the ion mobility on its size has also been explained from an entirely different point of view in terms of the dielectric friction.When an ion is placed in the polar medium, surrounding solvent is polarized due to the electric field of the ion.If the ion is displaced, the relaxation of the solvent polarization should take place to accommodate itself to the new position of the ion.The relaxation process is associated by an energy dissipation, and the corresponding drag force can be identified as the dielectric friction.The dielectric friction reduces with decreasing electrostatic field as the ion size increases.The two models, solvent-berg with the effective Stokes radii and the dielectric friction, explain at the same qualitative level the observed feature of the ion-size dependence of dynamics of small ions in polar liquids; the friction reduces with increasing ionic radii when the ions are sufficiently small, and then it begins to enhance as the radii further increase.As long as the general behaviour is concerned, all the theoretical treatments proposed so far, pure hydrodynamic [22][23][24], composite hydrodynamic and microscopic [25,26], could apparently have explained a qualitative nature of ion dynamics projected on the friction coefficient.However, the theories have left many important questions unsolved.Which of the above two models for ion dynamics is more realistic?If both processes are taking place, how are they coupled?What is the qualitative difference which distinguishes a solvent from other solvents with respect to ion dynamics?
Additional questions are raised when the moving ion is immersed in water instead of a simple dipolar liquid.The hydration phenomena observed by many different experiments show a remarkable variety in its dependence on size and sign of ions, which may not be characterized by a simple chemical model such as ionwater complex formation, or a solvent-berg model.Such a model seems to apply to very small monovalent ions like Li + and F − and multi-valent ions which make a stable hydration shell with a substantial lifetime.However, water molecules in the first hydration shell around ions with greater size are more mobile and disordered than those in bulk water.The behaviour which has been referred to as 'structure breaking' by Frank and Wen [27] and as 'negative hydration' by Samoilov [28] is due to the competition between two forces acting on water molecules in the shell, the isotropic electrostatic field from the ion and the hydrogen-bonding with other water molecules.Such characteristics of ion hydration is expected to produce non-linearity in the response of the electrostatic potential fluctuations by solvent to the ion field, which has strong dependence on size of ions [29,30].Recently, such hydration phenomena have been characterized theoretically based on the RISM theory in terms of the ion perturbation to the pair-correlation functions [31].The following questions may be naturally raised in conjunction of the ion dynamics in water.Does the dynamics of ions in water reflect such characteristics of the equilibrium structure of hydration?If it does, how and to what extent?Is there any qualitative difference in dynamics of ions in water compared to that in simple polar solvents, and what is the difference if any?In this section, we try to answer some of these questions.
Here, we address the problem from the response of collective density fluctuations in solvent, described in the previous section, to the ionic field [32,33].Since the Stokes and dielectric frictions originate basically from the energy dissipation associated with the translational and rotational motions of solvent, respectively, it is reasonable to ask how the ionic field couples to the collective density fluctuations, and/or how the two drag forces are related to the two collective modes.Since the translational and rotational motions of molecules are inherently coupled with each other in our atombased description of solvent dynamics, the theory is free from artifacts associated with the decoupling of those motions, which is inevitable in the earlier theories based on the explicit orientational coordinates.

Memory-function formulation of ion dynamics
Exploiting the standard fluctuation-dissipation theorem of Einstein, the friction coefficient ζ is related to the diffusion coefficient D as The diffusion coefficient in turn is related to the velocity auto-correlation function Z(t) by the Green-Kubo formula [34] The velocity auto-correlation function for an atomic ion in fluid is defined by where v u,z (t) denotes the z-components of the velocity of the ion at time t.(The subscript u signifies the solute ion.)The procedure of the projection-operator formulation leads to the memory-function equation for Z(t) [34]: where K(t) denotes the memory function.
Following the procedure proposed by Sjögren [35,36], we approximate the memory function as where K fast (t) denotes the rapidly decaying portion of the memory function due to the binary collisions, while K slow (t) represents the slow portion arising from correlated collisions.The fast portion can be well represented by the Gaussian ansatz, with the decay constant defined by In the expression, K(0) is the so-called "Einstein frequency", which can be calculated exactly from the information of the site-site intermolecular potential as well as of the density pair-correlation functions.K(0) is also an equilibrium quantity but with the three-particle correlation functions [32].
For the slow part of the memory function, we employ the mode-coupling theory, which leads to [32] where m is the mass of an ion, c uλ (k) denotes the site-site direct correlation function between the ion and solvent, and F u (k, t) represents the self-intermediate scattering function of the solute.We adopt the Gaussian approximation for F u (k, t), which is exact in the short-and long-time regime.F λµ (k, t) is the solvent site-site intermediate scattering function defined in the preceding section.f u (k, t) is an auxiliary function defined by denotes the intermediate scattering function of ideal gas.The velocity auto-correlation function Z(t) and its memory function K(t) can be obtained as a self-consistent solution of the above equations.

Velocity auto-correlation functions
The velocity auto-correlation functions (VACF) of alkali and halide ions in water obtained from the theory just described are plotted against time in figure 4. The general behaviour of VACF is high-lighted by the pronounced oscillation seen in the small ions, which disappears with increasing ion size.The oscillation is apparently  due to the vibrational motion of the ions.Similar pictures have been obtained by Rasaiah and co-workers based on the molecular dynamics simulation [37].This suggests the existence of a solvent cage around the ions, whose life-time is long enough to support the vibration.Remember, the ions which have pronounced oscillation in VACF are those classified as of 'positive' hydration by Samoilov [28].Therefore, the dynamic picture is well in harmony with the static one in this respect.The oscillation disappears in two physical causes; the size and mass.As the ion size increase, the ion-solvent interaction, whose essential nature is electrostatic, decreases, which loosen the solvent cage.Increasing mass makes the ion motion more 'inertial' with longer characteristic time.On the whole, the behaviour of VACF obtained from the theory is consistent with our intuitive picture for ion dynamics in water.

Diffusion constants
The diffusion constants calculated by equation (13) using VACF described above are depicted in figure 5 as functions of the ion radius, which is taken as half of the Lennard-Jones σ parameter.The behaviour is striking in a sense that it entirely breaks the Einstein-Stokes law, which predicts a monotonic decrease of the diffusion constant as the ion size increases.(Though we have not shown here, our theory in fact predicts the monotonic decrease, when the electrostatic interaction between the ion and solvent molecules is 'turned off.')Instead, our results exhibit just the opposite behaviour to the Einstein-Stokes law when the ion size is small, and turn into the same trend as the ion size is further increased.The physical cause of the anomalous behaviour seen in the diffusion constant we will discuss in detail later, because the diffusion constant contains essentially the same information as the friction coefficient.Here, we would like to make just two points.By comparing the plots for the cations (upper panel) and for the anions (lower panel), one can readily find that those curves do not fall on top of each other.The asymmetry with respect to the sign of charges comes from the asymmetry of charge distribution in a water molecule.As has been addressed in many ways, structure of water is characterized by the well-developed hydrogen-bond network, which has its origin in the

Figure 1 .
Figure 1.(a) Eigen frequencies as evaluated by diagonalizing ω 2k .Solid, dashed and dash-dotted lines give the eigen frequencies of the acoustic mode, OM-I and OM-II, respectively.(b) x O (solid line), x H 1 (dashed line) and x H 2 (dashdotted line) defined in the text corresponding to the acoustic mode.Notice that x H 1 = x H 2 holds in the whole kregion.(c) x O (solid line), x H 1 (dashed line) and x H 2 (dash-dotted line) corresponding to OM-I.Notice that x H 1 = x H 2 holds in the whole k-region.(d) x O (solid line), x H 1 (dashed line) and x H 2 (dash-dotted line) corresponding to OM-II.Notice that x H 1 = −x H 2 holds in the whole k-region.From (b) to (d), x O , x H 1 and x H 2 are normalized such that x 2 O + x 2 H 1 + x 2 H 2 = 1.

Figure 2 .Figure 3 .
Figure 2. The theoretical (solid lines) and MD simulation (circles) results of the longitudinal current spectra for the acoustic dynamics, C L,MM (k, ω).The results are plotted as a function of ω for k = n k min with k min = 0.3185 Å−1 and n = 1 ∼ 4, in arbitrary units.k min refers to the minimum accessible wave vector from MD simulation.

Figure 4 .
Figure 4.The normalized velocity autocorrelation functions of cations and anions in water.