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.

R ecent advances in nanoscale fabrication techniques have enabled the synthesis of nanofluidic systems with novel functionalities, 1−3 with applications to biotechnology, 4 filtration, 5−7 and computation. 8−10 For example, nanofluidics-based membranes have leveraged atomic level details like those of evolved biological membranes 11−18 to circumvent traditional trade-offs between permeability and selectivity that plague membrane technology. 19−22 While continuum-level hydrodynamic descriptions can remain accurate at scales of a few nanometers, enabling some general design principles to be deduced, 23−26 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. 27 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, 28,28−36 and the potential to incorporate selectivity for desalination or blue energy applications. 37−46 By computing the spatially resolved volumetric, charge, and species transport coefficients from equilibrium correlations 47−49 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, 50,51 as is consistent with observations that water density profiles are insensitive to wall fluctuations 52 and water−surface interactions are dominated by mean force fluctuations, 28 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 Å. 53,54 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 nm 2 . 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 w ≈ 1 g/cm 3 , computed via a temporal and spatial average over the central 2 nm of the channels. The channels were filled with N w = 1920 TIP4P/2005 water molecules with rigid geometries imposed using the SHAKE algorithm, 55 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. 52 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 dynamics where each particle i has mass m i , velocity v i , and experiences a friction ζ i , with forcing from interparticle interactions F r ( ) i N , and random noise R i . The random force is Gaussian with mean ⟨R i,α ⟩ = 0 and variance ⟨R i,α (t) R i′,α′ (t′)⟩ = 2k B Tζ i δ i,i′ δ α,α′ δ(t − t′) for each Cartesian coordinate α, where k B T 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, 57,58,58 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. 59 For all data presented here, we performed 5 independent simulations, each starting with an equilibration run for 5 ns with m i /ζ i = 2 ps, followed by a production run for 10−20 ns with m i /ζ 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. 60 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, i . We observe similar structures in both materials with interfacial layering of water that is consistent with previous simulations of neat water. 35, 58 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. 61−63 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, i , within 1.5 nm. We consider fluxes induced by a pressure differential, −ΔP x , 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 forces is the excess water flux, and z ( ) aB 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., jJ links the potential drop, −ΔΦ x , directly to the induced charge flux, j(z), while the off-diagonal elements are the so-called crossterms 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 as where particle i has velocity v i,x (t) and position z i (t) at time t and a static charge of c i , 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 and N b are respectively the average numbers of water molecules and all molecules in the bulk and A s is the surface area associated with the fluid−wall interface. The spatial dependence can be integrated out by defining total fluxes, such . 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. 47−49 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. 47 Mobilities are obtained as the long time slope of the correlation between time-integrated local and global fluxes with the correlation function 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. 28,35,64,65 In the presence of ions, the interfacial structure of water is altered and as a consequence the friction may change. In Figure 2 While the hydraulic conductivities deduced above are independent of a specific hydrodynamic model, they can be connected to continuum theory through the slip length l s . 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, v x = l s (∂v x /∂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 To apply a hydrodynamic interpretation, we consider only the region where a hydrodynamic description is expected to be  28 The inset illustrates the geometric interpretation of the slip length.

Nano Letters pubs.acs.org/NanoLett
Letter valid by defining the effective hydrodynamic interface as the location of the second water density peak in Figure 1(b). 26 where H hyd is the distance between hydrodynamic interfaces, and η is the estimated viscosity of the solution. This expression may be integrated to determine the hydraulic conductivity , and compares them against reported results for neat water. 28 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 l s has been observed in various contexts experimentally 32,66−69 and is understood to arise from a decoupling of structure and dynamics, though the mechanism is debated. 28,35,52,64,70,71 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. 72−74 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 jQ , 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 JQ 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 qJ being statistically zero at all points in the channel, consistent with . The pressure driven excess water mobility dQ , 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, dQ , with 0 dQ > corresponding to selective flux of water through the channel. We find a positive integrated value 0 dQ > 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 is where ρ 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 density This 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 . 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 l s ≫ H hyd /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 mobility The ratio of salt to total particle flux as a function of channel height is obtained as which is shown in Figure 4(b) normalized against the overall number fraction of ions in the bulk, f ( )/ salt K Cl tot = + + . 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 Nano Letters pubs.acs.org/NanoLett Letter 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, 75−77 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.
■ ASSOCIATED CONTENT