Phillip Helms,Anthony R. Poggioli,David T. Limmer
Intrinsic Interface Adsorption Drives Selectivity in Atomically Smooth Nanofluidic Channels
nanofluidics,desalination,membranes,interfaces,hydrodynamics
Specific molecular interactions underlie unexpected and useful phenomena in nanofluidic systems, but these require descriptions that go beyond traditional macroscopic hydrodynamics. In this letter, we demonstrate how equilibrium molecular dynamics simulations and linear response theory can be synthesized with hydrodynamics to provide a comprehensive characterization of nanofluidic transport. Specifically, we study the pressure driven flows of ionic solutions in nanochannels comprised of two-dimensional crystalline substrates made from graphite and hexagonal boron nitride. While simple hydrodynamic descriptions do not predict a streaming electrical current or salt selectivity in such simple systems, we observe that both arise due to the intrinsic molecular interactions that act to selectively adsorb ions to the interface in the absence of a net surface charge. Notably, this emergent selectivity indicates that these nanochannels can serve as desalination membranes.
Intrinsic Interface Adsorption Drives Selectivity in Atomically Smooth Nanofluidic Channels

Specific molecular interactions underlie unexpected and useful phenomena in nanofluidic systems, but these require descriptions that go beyond traditional macroscopic hydrodynamics. In this letter, we demonstrate how equilibrium molecular dynamics simulations and linear response theory can be synthesized with hydrodynamics to provide a comprehensive characterization of nanofluidic transport. Specifically, we study the pressure driven flows of ionic solutions in nanochannels comprised of two-dimensional crystalline substrates made from graphite and hexagonal boron nitride. While simple hydrodynamic descriptions do not predict a streaming electrical current or salt selectivity in such simple systems, we observe that both arise due to the intrinsic molecular interactions that act to selectively adsorb ions to the interface in the absence of a net surface charge. Notably, this emergent selectivity indicates that these nanochannels can serve as desalination membranes.

Recent advances in nanoscale fabrication techniques have enabled the synthesis of nanofluidic systems with novel functionalities, with applications to biotechnology, filtration, and computation. For example, nanofluidics-based membranes have leveraged atomic level details like those of evolved biological membranes to circumvent traditional trade-offs between permeability and selectivity that plague membrane technology. While continuum-level hydrodynamic descriptions can remain accurate at scales of a few nanometers, enabling some general design principles to be deduced, the continued development of nanofluidic devices is limited by a lack of understanding of emergent interfacial effects which are resolutely molecular in origin. With large surface to volume ratios, the properties of fluids confined to nanometer scales are determined in large part by a delicate interplay of interactions between the bounding surfaces and the working fluid. To understand and design nanofluidic devices, an approach that combines macroscopic and molecular perspectives is necessary. In this letter, we show how interfacial atomic structure affects the directed transport of an electrolyte solution in nanochannels made of atomically flat graphite (GR) and hexagonal boron nitride (BN) walls using molecular dynamics simulations unified with a contemporary perspective on hydrodynamics. These simple systems have been studied extensively because of their intriguing transport properties, such as anomalously high permeabilities in GR, and the potential to incorporate selectivity for desalination or blue energy applications. By computing the spatially resolved volumetric, charge, and species transport coefficients from equilibrium correlations we elucidate the importance of molecular interactions on nanofluidic device functionality. While from a continuum perspective, driving the solution with a pressure gradient should result in salt filtration or electric current only when the confining walls have a net charge, we discover that the intrinsic interfacial adsorption of ions can lead to streaming electrical currents and a novel, emergent desalination mechanism. We focus on the two systems illustrated in Figure 1(a), consisting of an aqueous solution of potassium chloride confined in nanochannels with walls of either BN or GR. Because hydrodynamics are insensitive to solid dynamics for stiff confining walls, as is consistent with observations that water density profiles are insensitive to wall fluctuations and water–surface interactions are dominated by mean force fluctuations, the walls are held fixed. Because of the experimental similarity between the structure of BN and GR lattices, we spaced atoms and lattice layers identically, with interatomic and interlayer spacings of 1.42 and 3.38 Å. Each wall has three layers, using AA′ and AB stacking for BN and GR, respectively, to match their equilibrium structures, with lattice unit cells repeated 8 and 13 times in the x and y directions for a cross-sectional surface area of nearly 9 nm2. The walls were separated such that the spacing between the center of mass of the innermost wall layers was H ≈ 5.7 nm, with the channel width adjusted to ensure a bulk water density of ≈ 1 g/cm3, computed via a temporal and spatial average over the central 2 nm of the channels. The channels were filled with Nw = 1920 TIP4P/2005 water molecules with rigid geometries imposed using the SHAKE algorithm, potassium ions, and chloride ions, resulting in a nearly 1 M electrolyte solution. Previous work evaluating the friction of water on both BN and GR found the predominate contribution was thermodynamic in origin, characterized by the mean force fluctuations between the liquid and solid (rather then how those force correlations decay dynamically). Recent work has found that the water density profiles, and thus mean force fluctuations, are indistinguishable if the surfaces are allowed to fluctuate or not. This is likely because both BN and GR are stiff lattices and, as a consequence, we do not expect the freezing of the lattices to impact our observations. We evolved this system according to underdamped Langevin dynamicswhere each particle i has mass mi, velocity vi, and experiences a friction ζi, with forcing from interparticle interactions , and random noise Ri. The random force is Gaussian with mean ⟨Ri,α⟩ = 0 and variance ⟨Ri,α(t) Ri′,α′(t′)⟩ = 2kBTζiδi,i′δα,α′ δ(t – t′) for each Cartesian coordinate α, where kBT is Boltzmann’s constant times temperature. Periodic boundary conditions were imposed in all three spatial dimensions, with a vacuum layer in the z direction of 5 nm to ensure no interaction between periodic images of the channel. Intermolecular Lennard-Jones forces were chosen from literature-reported values to reproduce the solubility of ions in water and match the ab initio equilibrium fluid structure in BN and GR nanochannels, with Lorentz–Berthelot mixing rules defining heteroatomic interactions. Additionally, water molecules, charged ions, and the BN wall atoms interacted with Coulomb potentials, where boron and nitrogen atoms have charges of ±1.05e, with e being the elementary charge, using an Ewald summation as implemented in LAMMPS. For all data presented here, we performed 5 independent simulations, each starting with an equilibration run for 5 ns with mi/ζi = 2 ps, followed by a production run for 10–20 ns with mi/ζi = 10 ns at a temperature of 298 K. In all plots, lines represent averages and shaded bands represent the standard deviation for the 5 simulations. All scripts used to produce these results and the raw data are openly available. Figure 1(b) shows the equilibrium particle number densities, ρi(z) for water, potassium, and chloride, with i = {w, K+, Cl–}, in the BN and GR channels, relative to their bulk values, . We observe similar structures in both materials with interfacial layering of water that is consistent with previous simulations of neat water. The distribution of ions near such interfaces is known to be highly dependent on ion species, and the profiles shown are consistent with previous simulations. A dense layer of pure water accumulates near the wall, with the molecules oriented such that they induce a small local negative charge. The next layers are enriched in alternating concentrations of potassium and chloride ions, with depletion (accumulation) of water molecules accompanying potassium (chloride) enrichment. The two materials differ slightly, with a higher water density in the first layer of BN resulting in layering with higher amplitude in BN compared to GR, though in both systems the layering in the density decays to its bulk value for each species, , within 1.5 nm. We consider fluxes induced by a pressure differential, −ΔPx, imposed electrostatic potential drop, −ΔΦx, or water chemical potential differential, −Δμx, with subscripts denoting application in the x direction parallel to the walls, and limit ourselves to small driving strengths. In this limit, linear response theory dictates that induced local fluxes are linearly dependent on driving forceswhere q(z) is the volumetric flow, j(z) is the charge flux, d(z) is the excess water flux, and are the spatially dependent mobilities. The excess water flux d(z) represents the local water flux relative to what would be predicted from the bulk water density and the local total flux of water and ions, and it is considered here because it is particularly relevant for desalination. The diagonal elements of the mobility matrix link a given forcing directly to its conjugate flux, e.g., links the potential drop, −ΔΦx, directly to the induced charge flux, j(z), while the off-diagonal elements are the so-called cross-terms linking, for example, an induced charge flux to an applied pressure differential. The total fluxes include the total volumetric flow Q, charge flux J, and excess water flux D. We index mobilities by the local induced flux a and total flux B directly conjugate to a particular forcing. The local fluxes are defined microscopically aswhere particle i has velocity vi,x(t) and position zi(t) at time t and a static charge of ci, and δi,w is a Kronecker delta that returns 1 if particle i is a water molecule and 0 otherwise. The bulk water mole fraction is defined as , where and Nb are respectively the average numbers of water molecules and all molecules in the bulk and As is the surface area associated with the fluid–wall interface. The spatial dependence can be integrated out by defining total fluxes, such as , with analogous definitions for J and D. Total channel conductivities can be evaluated as , resulting in total flux linear response relations such as . While the integrated conductivities must obey Onsager reciprocal relations, , mobilities are under no such constraint. It is possible for . Rather than attempting to calculate mobilities directly via nonequilibrium simulations, we use fluctuation–dissipation relations in order to obtain transport coefficients from equilibrium flux correlations. This allows us to avoid running separate nonequilibrium simulations for each term in the mobility matrix, and ensures the validity of linear response. We adopt the Einstein–Helfand approach over the Green–Kubo method, as recent work has demonstrated its enhanced statistical efficiency. Mobilities are obtained as the long time slope of the correlation between time-integrated local and global fluxeswith the correlation functionvolume V = AsH, and brackets representing an equilibrium average. Similarly, conductivities can be obtained using correlations between global fluxes, with . Previous work has demonstrated that while equilibrium structures suggest minor differences between water in BN and GR nanochannels, the fluid dynamics are strikingly different. This results in significant differences in friction at the fluid-wall interface and thus channel permeabilities. In the presence of ions, the interfacial structure of water is altered and as a consequence the friction may change. In Figure 2(a), we show the integrated global flux correlation function KQQ as a function of time for both nanochannels. After approximately 200 ps, the correlation functions approach a linear dependence on time and their slopes give the hydraulic conductivities as = 18.0 ± 9.2 mol nm5 kJ–1 ns–1 and = 106 ± 40 mol nm5 kJ–1 ns–1, which differ by nearly an order of magnitude. While the hydraulic conductivities deduced above are independent of a specific hydrodynamic model, they can be connected to continuum theory through the slip length ls. In contrast to the no-slip condition typically applied in macroscopic contexts, which specifies that the fluid velocity vanishes at the walls, the confinement and enhanced interfacial importance in nanofluidic applications typically require the finite-slip condition. This condition specifies that the velocity at the wall is proportional to the shear strain at the wall, vx = ls(∂vx/∂z)|z=0. The slip length is interpreted geometrically as the distance beyond the interface where the extrapolated flow profile is zero, illustrated in Figure 2(b). To apply a hydrodynamic interpretation, we consider only the region where a hydrodynamic description is expected to be valid by defining the effective hydrodynamic interface as the location of the second water density peak in Figure 1(b). At this distance, microscopic density correlations have decayed and the fluid is well described as a continuous medium. The Poiseuille solution for the hydraulic mobility in the presence of a finite slip length is given bywhere Hhyd is the distance between hydrodynamic interfaces, and η is the estimated viscosity of the solution. This expression may be integrated to determine the hydraulic conductivitywhich allows us to relate the measured values of in GR and BN to the corresponding slip lengths provided η is known. Here, we use a viscosity of η = 1.0 mPa s, obtained by interpolating literature values.Figure 2(b) indicates the resulting slip lengths, and , and compares them against reported results for neat water. With the slip being approximately an order of magnitude larger in GR than in BN nanochannels, it is clear that the qualitative results do not change significantly with the addition of salt. The material-dependency of ls has been observed in various contexts experimentally and is understood to arise from a decoupling of structure and dynamics, though the mechanism is debated. Quantitatively, our simulations also suggest a decrease in slip as salt is added, consistent with observations for slip on hydrophobic surfaces, where increasing fluid–wall friction results as a consequence of enhanced equilibrium force fluctuations from the heterogeneous solution. Detailed insight into the differences in transport characteristics between these nanochannels can be obtained by computing the spatially dependent hydraulic mobility using eq 4. The results of this calculation are shown in Figure 3(a). We also show the corresponding profiles calculated from eq 6 for comparison to the macroscopic theory. As expected for the conductivity, we observe approximately an order of magnitude difference between the peaks in the hydraulic mobilities in the BN and GR nanochannels. The mobility profile is nearly flat for GR and exhibits a slight curvature for BN, indicative of the differences in slip. In the boundary region, the mobility profile qualitatively mimics the fluid density profile with greater (lesser) flux coinciding with density peaks (troughs). We find that the molecular interfacial structure also affects the cross-terms in the mobility matrix in eq 2. The streaming mobility , which quantifies the electrical current profile produced by applying a pressure differential, is shown in Figure 3(b) for both systems. We observe the emergence of three layers of electrical current of alternating sign near the fluid-wall boundary, and no net current in the bulk. Because the applied pressure produces particle flux in all regions, the alternating current is caused by ion density localization at the interface, with positive (negative) current where potassium (chloride) ions are enriched. These interfacial effects decay away from the wall more slowly than those observed with the hydraulic mobility, with net charge flux penetrating into the hydrodynamic region defined by the hydraulic mobility. By integrating the mobility across the channel, we find that the streaming conductivity is indistinguishable from zero for both materials, indicating no net ionic transport. Though not shown, our calculations verify the lack of symmetry between cross-term mobilities, with being statistically zero at all points in the channel, consistent with while maintaining . The pressure driven excess water mobility , is shown in Figure 3(c) as computed using eq 4. This quantity is related to the desalination capabilities of a nanochannel, with its magnitude determined by the channel’s selectivity and permeability. This transport is summarized by the integrated mobility, , with corresponding to selective flux of water through the channel. We find a positive integrated value for both materials, demonstrating water selectivity, corresponding to salt rejection of approximately 25%, to be illustrated later. The spatial dependence of the cross-term mobility profiles can be understood via a combination of microscopic and macroscopic perspectives. The streaming mobility may be evaluated microscopically as a product of density profiles and the hydraulic mobility. For the streaming mobility this iswhere ρtot(z) = ρw(z) + ρK+(z) + ρCl–(z). Though a common decomposition in macroscopic hydrodynamics, this is a nontrivial statement when considering the microscopic mobilities. The red dashed line in Figure 3(b) shows this estimate agrees well with estimate using eq 4. The same functional decomposition holds for excess water flux, obtained from the product of the hydraulic mobility and excess water densityThis decomposition is shown in the red dashed line in Figure 3(c). Both of these decompositions follow directly from the Langevin equations of motion. While the excess water mobilities for both materials are qualitatively similar because of their similar density and hydraulic mobility profiles, the quantitative difference arises due to the differences in magnitude of the hydraulic conductivity. The first contact layer is nearly salt-free, so while interfacial friction slows pressure driven transport, the high water purity gives a large peak in excess water mobility. There is a second excess water mobility peak near the second water density peak. The enrichment and depletion of chloride and potassium, respectively, brings the overall salt density close to its bulk value and leaves an excess concentration of water where the hydraulic mobility also peaks. The molecular dynamics calculations suggest that the transport properties of the nanochannel can be decomposed as a sum of a molecular interfacial component, and a continuum bulk component. The interfacial component depends sensitively on specific molecular interactions as they manifest in nonuniform density profiles. Beyond the domain of those density correlations, which for these channels extend around 2 nm into the channel, the transport is well described by Poiseuille flow with a large slip length. This decomposition allows us to infer the height dependence of the channel’s selectivity and permeability. We can calculate the size dependent conductivity using an integrated mobility , where we employ inversion symmetry to integrate over only half of the channel. These conductivities are shown for in Figure 4(a) normalized against . The red regions in Figure 4 indicate system sizes which would lead to overlapping interfacial regions, for which our decomposition is likely invalidated. This boundary is larger than the hydrodynamic region shown in Figure 2 because interfacial effects on ion density extend into the expected hydrodynamic region. Because the hydraulic mobility profile is nearly flat in the hydrodynamic region, which is expected when ls ≫ Hhyd/6, the permeability increases linearly with channel height, which is slower than anticipated from traditional hydrodynamics with no slip. A similar approach can be used to compute the dependency of the water selectivity on the height of the channel. To compute the selectivity, we first can determine a pressure driven salt mobilityThe ratio of salt to total particle flux as a function of channel height is obtained aswhich is shown in Figure 4(b) normalized against the overall number fraction of ions in the bulk, . This provides a direct measurement of the size dependence of the nanochannel selectivity. Consistent with the inference from the excess water mobility, the salt flux is suppressed relative to its expected value from the bulk concentration of ions and the total channel conductivity. We find that BN and GR nanochannels have effectively identical selectivities, primarily because of their similar equilibrium fluid density distributions and qualitatively similar hydraulic mobility profiles. For the nanochannel size and ion concentrations considered here, the flux of salt ions is reduced by approximately 25%, while shrinking the nanochannel until interfacial regions overlap at around 2 nm could provide a reduction of around 50%. Due to the intrinsic interfacial absorption of ions to the interface and their resultant suppressed mobility, as the nanochannel size is decreased its selectivity is enhanced. An optimal desalination device must separate ions from water with both high selectivity as well as high permeability, and these phenomenological channel scaling observations suggests that for both BN and GR this optimum is between 2 and 5 nm. This mechanism of selective transport, and the ability of the channel to separate salt from water, is a result of an interplay between local molecular interactions that drive ions to the fluid–solid boundary in the absence of a net surface charge of the substrate. These molecular interfacial features established a nonuniform fluid composition across the channel that, when combined with a spatially resolved evaluation of the hydraulic mobilities, provide a complete description of the transport within the nanochannel. The promise of this mechanism for desalination technology is strikingly enhanced when this water selectivity is coupled with the anomalously high permeability of GR nanochannels. This framework is general and can be used to understand and engineer other functionality in nanofluidic systems. Employing recent generalizations of response theory, our approach could be extended outside the regime of linear response to provide insight into performance at high driving strengths and between multiple driving forces.