Current fluctuations in nanopores: The effects of electrostatic and hydrodynamic interactions

Using nonequilibrium Langevin dynamics simulations of an electrolyte with explicit solvent particles, we investigate the effect of hydrodynamic interactions on the power spectrum of ionic nanopore currents. At low frequency, we find a power-law dependence of the power spectral density on the frequency, with an exponent depending on the ion density. Surprisingly, however, the exponent is not affected by the presence of the neutral solvent particles. We conclude that hydrodynamic interactions do not affect the shape of the power spectrum in the frequency range studied.

Hydrodynamic interactions have a strong influence on the dynamics of Brownian particles suspended in a solvent, producing self-organized states, nonlinear dynamics, and synchronization [1][2][3][4]. Hydrodynamic interactions between objects decay slowly. Similar to the electrostatic potential, the strength of the hydrodynamic interactions in bulk is inversely proportional to the distance between the particles [5]. Moreover, the effects of hydrodynamic interactions are extremely sensitive to geometric confinement. Density perturbations in a fluid between stationary confining walls give rise to a long-time tail in the velocity autocorrelation function of colloidal particles [6]. Experiments show that hydrodynamic interactions even become independent of the distance between particles inside small pores [7]. These hydrodynamic effects have a pronounced effect on the dynamics of larger molecules, such as DNA, translocating through a nanopore [8,9]. Whereas the effect of hydrodynamic interactions on colloidal particles and polymer dynamics has attracted a lot of attention over the past decades, the effect of hydrodynamic interactions on ion dynamics remains largely unexplored.
The combination of experimental measurement and molecular modeling of the power spectral density constitutes a promising technique to study ion motion in unprecedented detail [10]. For example, the power spectrum can be used to study the microscopic properties of nanofluidic systems, such as the adsorption of molecules on the walls of a nanometer-scale cavity [11]. Recently, we showed that ion correlations at high particle density produce a power-law spectrum at low frequency, with an a e-mail: ramin.golestanian@physics.ox.ac.uk b e-mail: douwe.bonthuis@physics.ox.ac.uk exponent depending on the ion density [12]. Experiments show that the power spectrum of the ionic nanopore current S(ω), with ω = 2πf being the frequency, typically follows a power law S(ω) ∝ 1/ω α with α ≈ 1, which is referred to as pink noise, or 1/f noise [13][14][15]. The appearance of pink noise in nanopore current measurements is ubiquitous; it is found in a variety of systems, from protein channels and flexible synthetic pores [16][17][18], to solid-state conical pores [19]. The molecular origin of the low-frequency pink noise has been debated for decades [20][21][22]. However, theoretical analysis of the power spectrum including the multi-body interactions between the ions and the effect of hydrodynamics remains challenging. Therefore, although hydrodynamic interactions are usually present in experimental studies, their effects on the frequency dependence of the power spectral density are unknown. The situation has changed with the recent advance of fast and versatile molecular simulation techniques, which now allow a systematic computational investigation.
In this manuscript, we present a Langevin dynamics simulation study of the ionic current through a nanometer-scale pore filled with an electrolyte, using the Espresso molecular dynamics package [23]. The electrolyte is modeled by ions in an explicit solvent. For the solvent, we use a coarse-grained description of neutral, nonpolar Lennard-Jones particles. To systematically study the effect of hydrodynamic interactions, we vary the density of both the ions and the solvent particles independently. We calculate the power spectral density of the ion current and compare the results with simulations without solvent and with a linearized mean-field theory of ion currents without hydrodynamic interactions. Whereas an increase in the ion density directly causes a power-law behavior of the power spectrum at low frequency, introducing hydrodynamic interactions by increasing the solvent density does not have the same effect.

Simulation model
We use the simulation package Espresso [23] to set up Langevin dynamics simulations of a nanopore filled with a mixture of monovalent positive and negative ions and neutral solvent particles (Fig. 1). The Langevin equation for particle i is expressed as where ξ i (t) is the stochastic force satisfying ξ i (t) · ξ j (t ) = 6k B T γδ ij δ(t − t ), u i and m i = 1 k B T τ 2 /Å 2 denote the velocity and the mass, respectively, and F i is an external force applied to the particle. The thermal energy equals k B T and we use γ = 1 k B T τ/Å 2 . By using an equal and arbitrary mass for all particles, m i is incorporated in the time scale τ . Interactions between pairs of particles are modeled by the sum of a Coulomb and a Weeks-Chandler-Andersen (shifted Lennard-Jones) potential V ij (r ij ), where r ij denotes the distance between particles i and j, Q i is the charge of particle i in units of the elementary charge e and l B = e 2 /(4πεε 0 k B T ) is the Bjerrum length. The Lennard-Jones radius is given by σ ij and ij represents the interaction strength.
The Lennard-Jones interaction is truncated at r ij = 2
The electric field on the ions is represented by a force applied inside the pore in the x direction, which is the direction along the length L of the pore (Fig. 1), The electric field is varied between E = 0.3 and 1.6 k B T /(eÅ). The simulations are performed in a cylindrical nanopore with a radius ranging from 19Å to 30Å, permeating a rigid membrane with a width of W = 96Å and a length L = 48Å (Fig. 1). We use an increasing solvent density of C s = {5.5 · 10 −4 , 2.1 · 10 −3 , 5.5 · 10 −3 }Å −3 , each of which is well below the bulk freezing density of the Weeks-Chandler-Andersen fluid model. Compared to the molecular density of water, the maximum density corresponds to a coarse-grained force field where each particle represents approximately 6 water molecules. A smooth surface induces a crystalline order in the fluid, over a range depending on the molecular properties of the liquid [24]. For a fluid consisting of identical Lennard-Jones spheres, the induced order propagates over a distance larger than our simulation box. Therefore, to prevent crystallization of the Lennard-Jones fluid, we perturb the uniform membrane surface by randomly removing half of the particles from the outer layer. The membrane particles are frozen, and for the membrane-ion, membrane-solvent, ion-solvent, ion-ion and solvent-solvent interactions we use ij = 2 k B T and σ ij = 4.7Å. For the membrane and solvent particles we use Q i = 0, and for the ions we use Q i = ±1. The ion concentration is varied according to When the motion of particle i perturbs the surrounding solvent, the hydrodynamic signal diffuses at a rate governed by the kinematic viscosity ν. For hydrodynamic interactions to occur, this viscous momentum must diffuse much faster than the particle itself. The relation is governed by the Schmidt number Sc = ν/D, with D being the diffusion coefficient of the solvent particles. To verify that the coarse-grained solvent particles produce hydrodynamic interactions in the strongly confined environment of the nanopore, we simulate a pressure difference across the length of the with R being the radius of the pore, F C s = ∇ p being the pressure gradient across the length of the pore, which in our case is derived from the uniform applied force F = 0.8 k B T /Å on the solvent particles inside the pore, which have mass m i = 1 k B T τ 2 /Å 2 and number density C s . We show the velocity profile in Fig. 2(a) for the lowest solvent density C s = 5.5 · 10 −4Å−3 . The fit of Eq. (4) yields ν = 900Å 2 /τ , which in combination with D = 1Å 2 /τ yields Sc = 900. As the Schmidt number at higher solvent densities is even higher, all our simulations satisfy the condition for hydrodynamic interactions. Apart from propagation by viscous momentum diffusion, hydrodynamic interactions are transmitted by sound wave propagation. In an incompressible fluid, the sound velocity is infinite, and the viscous momentum diffusion is solely responsible for the time evolution of the hydrodynamic interactions. The compressibility of our model solvent is finite, however, depending on the solvent density C s , which might have implications for the hydrodynamic interactions [25]. We calculate the isothermal compressibility from the pressure p as a function of solvent density C s in separate bulk simulations using κ T = (d ln(C s )/dp) T , with p being the pressure, see Fig. 2(b). The compressibility is varied over two orders of magnitude as we change the solvent density. Nevertheless, the compressibility of water, equal to κ T = 2Å 3 /(k B T ), is still a factor 5 below the compressibility of our highest-density solution. We quantify the effect of the compressibility by calculating the sound velocity u s = γ/ (m i C s κ T ), with γ being the heat capacity ratio, which is of order γ ∼ 1 in a liquid. The sound velocity increases drastically when we change the solvent density in our simulations, from u s = 3Å/τ at the lowest density to u s = 239Å/τ at the highest density. The importance of the compressibility effects is estimated from the Mach number Ma = k B T /m i /u s , where the estimate of the thermal velocity k B T /m i is being used as the typical velocity of the particles. As Ma is well below 1 in all our simulations, the compressibility is not expected to have a large effect on the hydrodynamic properties [25].
The stochastic force ξ in the Langevin dynamics simulations provides a truncation length beyond which ξ exceeds the force due to hydrodynamic interactions [25]. Quantification is complicated, because the truncation length depends on the magnitude of the force from which the hydrodynamic interactions originate. As the interparticle forces in the system reach very high values, however, a part of the long-ranged hydrodynamic interactions will be preserved.

Linearized mean-field theory
We derive a theoretical description of the noise spectrum of the ionic current following our previous analysis [12]. The expression for the power spectral density S(ω) is derived for monovalent ions in implicit water. Therefore, the following derivation does not include the effect of hydrodynamic interactions. Comparison with the simulation results allows us to study the effect that the hydrodynamic interactions in the simulations have on the power spectral density. Ion-ion correlations, which are responsible for the low-frequency power-law increase of the power spectrum at high ion density [12], are also absent from this theoretical model. We consider a system consisting of a cylindrical nanopore of length L and radius R connecting two reservoirs (Fig. 1), and calculate the flux density J ± (x, t) of positive and negative ions inside the nanopore, with x denoting the position in three dimensions and t denoting the time. The ion concentrations C ± (x, t) are governed by the continuity equation, The corresponding flux densities J ± (x, t) are given by the Nernst-Planck equation, where E (x, t) is the applied electric field, e denotes the elementary charge, and η ± (x, t) denotes the thermal noise that accounts for fluctuations in the environment; most importantly the effect of the implicit water on the ion dynamics. From here, we switch to index notation where α, β, and γ correspond to the three components of our coordinate system. To simplify the notation, we assume D + = D − = D and η + x, t = η − x, t = ηx, t. After applying a standard Fourier transform to Eqs. (5) and (6), we find with ... denoting the Fourier transform, q being the wave vector, and ω being the frequency. Rewriting Eq. (7) leads to where M αβ q, ω denote the matrix components 1588 The European Physical Journal Special Topics Combining Eqs. (8) and (9) and solving for J ± γ , we find with det(M ) denoting the determinant of M . Within the geometry of the pore, there is one parallel ( ) direction, and two equivalent perpendicular (⊥1, ⊥2) directions, see Fig. 1(a). The electric field is nonzero only in parallel direction E = (0, 0, E ). Therefore, the flux in the parallel direction becomes with q ⊥1 and q ⊥2 being the two independent wave vectors in the plane of the membrane. As the random force is applied to every individual particle, the power spectrum of the thermal noise in our implicit-solvent model is proportional to the ion concentration inside the pore, with C V = N /(πR 2 L) being the average number of ions per unit volume in the pore, which is proportional to the bulk ion concentration C i , but depends nontrivially on the radius R, the length L, the electric field, and the interionic interaction potential. Introducing short-hand notation, we derive from Eqs. (11)-(12) with q 2 = q 2 ⊥1 + q 2 ⊥2 + q 2 . The two-sided power spectral density S (ω) of the current I (t) defined on the domain 0 < t < T is given by the limit of T → ∞ of which can be written as We rewrite I (t) as the integral of the current density J + (x, t) − J − (x, t) at a given position in the direction of x over the lateral surface area A of the pore, Some mathematical manipulation yields in the limit T → ∞, with Λ being the small-scale cut-off length, introduced because of the finite particle size. The Fourier-transformed area function in Eq. (17) is given by The integral in Eq. (18) is performed over the lateral area A of the pore, which is approximately circular. However, because our cylindrical direct space does not map exactly to a cylindrical reciprocal space, we use two different approximations to calculate the integral (Fig. 1(c)). First, integrating over a square of sides 2R gives Alternatively, integrating over a circle of radius R gives with x ⊥ = x 2 ⊥1 + x 2 ⊥2 and θ = arctan(x ⊥2 /x ⊥1 ) being the cylindrical coordinates, q ⊥ = q 2 ⊥1 + q 2 ⊥2 and φ = arctan(q ⊥2 /q ⊥1 ) being the polar coordinates in reciprocal space, and J 1 being the first order Bessel function of the first kind. The primary difference between Eqs. (19) and (20) is the amplitude of the calculated noise spectrum [12]. Contrary to the circular area, however, the square area can be mapped directly to reciprocal space, enabling a straightforward evaluation of Eq. (17). Therefore, we use Eq. (19) for all the curves in the present paper. Together with Eqs. (13) and (19), Eq. (17) is solved numerically to get the linearized mean-field prediction of S(ω).

Results and discussion
We calculate the power spectral density of the ion current in simulations with three different solvent densities C s (Fig. 3). At low solvent density, the curves exhibit a transition around ω = 0.1/τ , similar to the implicit solvent case [12]. With increasing solvent density, the transition becomes less pronounced due to an increase in the highfrequency noise level. Surprisingly, the increasing solvent density does not induce any alteration of the power spectral density at low frequency, even for a tenfold increase in solvent density. We verify the effect of increasing ion concentration on the power spectral density in the presence of hydrodynamic interactions (Fig. 4). The amplitude of the noise increases with increasing ion concentration, and the transition frequency shifts to slightly lower values. Most strikingly, however, is the change of the behavior at low frequency. The power spectral density exhibits a power law, with an exponent that increases sharply with increasing ion concentration. These results are similar to the results found in simulations with implicit solvent [12]. However, as the curves in Fig. 4 extend to higher ion concentrations than those treated previously, the new results show that the increase in the exponent of the power law continues, reaching a = 0.4 at an ion concentration of C i = 2 · 10 −3Å−3 .
To test the effect of the hydrodynamic interactions, we fit the linearized mean field theory, which does not take hydrodynamic interactions into account, to the curves in Fig. 4. Apart from the low-frequency power-law dependence, which is caused by ionion correlations [12], the simulated curves are well described by the implicit-solvent model. Remarkably, it is not necessary to take hydrodynamic interactions into account to describe the power spectrum of the ionic current through an electrolyte-filled pore.
At low frequency, we fit the exponent a of the power law S(ω) ∼ ω −a for the curves shown in Figs. 3 and 4. We fit the noise spectra for log 10 ω < −1.8 and discard the lowest frequency data points because of their statistical uncertainty. The exponent is shown in Fig. 5 as a function of ion concentration C i for fixed solvent show the simulation results and the shaded region represents the standard deviation that we get when applying a block-averaging method to improve the readability. The dashed lines represent the fits derived from the linearized mean field theory (Eqs. (13), (17) and (19)), with parameters taken from the simulations: applied electric field E = 1.6 kBT /(eÅ), pore radius R = 25Å, diffusion coefficient D = 1Å 2 /τ , and the small-scale cutoff length is set to a value of the order of the ion size, Λ = 2.5Å, for all curves. The solid black lines indicate the fit with S ∼ 1/ω a .    concentration (top panel) and as a function of solvent concentration C s at fixed ion concentration (bottom panel). Whereas the exponent increases sharply as a function of the ion density, increasing the solvent concentration has no effect. Because the charge is the only difference between an ion and a solvent particle, we conclude that electrostatic interactions cause the increasing exponent. Hydrodynamic interactions, despite having a similar long-ranged spatial dependence, do not have the same effect.
We perform an extra simulation without solvent particles (C s = 0), and compare the power spectra of simulations with and without explicit solvent directly in Fig. 6. Clearly, the curves have the same frequency dependence over the entire frequency range, confirming the results of the preceding paragraphs.
Finally, we study the dependence of the power spectral density on the pore radius R and the applied electric field E . In Fig. 7, we show that the linearized mean-field theory -derived for implicit solvent -captures the dependence on the pore radius and the electric field without further fit parameters, for all values of R and E studied.

Summary and conclusions
We present a systematic numerical investigation of the effects of hydrodynamic and electrostatic interactions on the power spectral density of ionic currents in nanopores using an explicit coarse-grained solvent. We find that an increase in ion concentration at fixed solvent density leads to a power-law behavior at low frequency with an exponent increasing with ion density. The power-law frequency-dependence of the power spectrum is in line with our previous findings in simulations with implicit water, where the nonzero exponent was shown to be caused by ion-ion correlations. The exponent reaches a = 0.4 at an ion concentration of C i = 2 · 10 −3Å−3 . Hydrodynamic interactions influence the power spectral density at high frequency. In particular, the transition in the power spectral density becomes less pronounced with increasing solvent density. At low frequency, however, the hydrodynamic interactions have no effect, which is surprising in the view of the large influence of hydrodynamic interactions on the dynamics of colloids and polymers under confinement. Note, however, that the solvent used in the present study has a higher compressibility than water, and that the Langevin noise provides a truncation distance, which might influence the hydrodynamic interactions. The linearized mean-field theory without hydrodynamics, which has been derived in our previous work [12], can be used to describe simulation results with hydrodynamic interactions equally well. Instead, inclusion of electrostatic ion-ion correlations is paramount to describe the low-frequency power-law behavior as a function of ion density. Although a direct comparison with experimental results is not yet feasible, we show that the combination of simulations and analytical work provides a promising framework for the systematic investigation of experimental noise spectra.