Unexpected behaviors in molecular transport through size-controlled nanochannels down to the ultra-nanoscale

Ionic transport through nanofluidic systems is a problem of fundamental interest in transport physics and has broad relevance in desalination, fuel cells, batteries, filtration, and drug delivery. When the dimension of the fluidic system approaches the size of molecules in solution, fluid properties are not homogeneous and a departure in behavior is observed with respect to continuum-based theories. Here we present a systematic study of the transport of charged and neutral small molecules in an ideal nanofluidic platform with precise channels from the sub-microscale to the ultra-nanoscale (<5 nm). Surprisingly, we find that diffusive transport of nano-confined neutral molecules matches that of charged molecules, as though the former carry an effective charge. Further, approaching the ultra-nanoscale molecular diffusivities suddenly drop by up to an order of magnitude for all molecules, irrespective of their electric charge. New theoretical investigations will be required to shed light onto these intriguing results.

(2) including only hard-sphere (HS) and hydrodynamics effects (solid blue lines), only exclusion or enrichment effects (green solid lines), and the full Eq. (2) (solid red line). The green line is computed from Eq. (2) with rs = 0, corresponding to neglecting HS and hydrodynamic effects. A much better fit can be achieved tuning both the molecule charge, expressed in the parenthesis, and volume [red dash line]. Panels (e, f ) show diffusivity data for two neutral species. As in the previous panels, blue solid lines are computed from Eq. (2) with q = 0, but a much better fit can be obtained considering the molecules charged. Each of the data points denotes the average of three individual replicates, error bars are ± sd.

Supplementary Note 1: Geometry and Fabrication Process
The membranes were fabricated using standardized industrial processes at a commercial manufacturer as previously described [8]. Each silicon chip had dimensions of 6 by 6 by 0.750 mm (length, width, and thickness) and enclosed large numbers of monodispersed nanochannels (over 50,000 mm −2  The capping dielectric stack, and hence the internal stress associated with these dielectrics, was consistent for nanochannel sizes between 2.5 and 40 nm. The stack thickness was far greater than the nanochannel height, leading to a nearly-identical structure and geometry for all devices. The thickness was changed for the 250 nm channels, as they required a different integration scheme. The stress was high enough to keep the ceiling of the channels tensile in all cases. Following the sacrificial etching process, isopropyl alcohol was substituted for water prior to drying to minimize surface tension and prevent channel collapse. The nanochannel walls were asymmetrically and negatively charged, as the potentials of -28.4 and -18.5 mV were indirectly measured with a zeta potential analyzer (Delsa Nano C, Beckman Coulter, Inc., CA) for the silicon and silicon nitride surfaces, respectively. The resulting slit-channels, presenting a highly defined and precise geometry with sub-nanometer size tolerances, were parallel to the membrane surface and orthogonal to the inlet and outlet microchannels; a design promoting both high channel density and physical robustness (See Figure 1a) [6]. this method achieved sufficient nanochannel density to permit significant outflow at a rate easily and reliably detectable by common instrumentation [8], overcoming a primary limitation inherent to similar nanoscale systems [3,7].

Gas-Flow Analysis
Silicon membranes were tested by measuring the gas-flow of nitrogen gas (Matheson Tri-Gas, NJ) though the nanochannels. This non-destructive, fast, and highly repeatable method was developed to ensure the fabrication quality of the whole population of silicon channels, addressing a drawback of SEM and TEM imaging inspection (limited to a small portion of one nanochannel). The nitrogen gas was pressurized at 15 psi by an automated pressure controller (PC3, Alicat Scientific, AZ) and forced to flow though the membrane, which was clamped by a custom-made, stainless steel holder.
Leakage within the membrane holder was prevented by a series of silicon o-rings (Apple Rubber, CA). The gas flux was measured by three different mass flow meters (M500SCCM, M10SCCM, M05SCCM, Alicat Scientific, AZ), whose employment depended on the sensitivity required (by dimension of the nanochannels). Over 30 measurements were taken once the pressure and flux reached a steady state for each membrane and averaged. To ensure the quality of the reading and the absence of leakage in the system, a blank test was performed every 15 measurements. This blank test employed a silicon membrane with the same dimension and structure as the membrane previously described, but lacking channels, allowing detection of any gas leakage from non-perfect seals. For calibration and data quality, any measurements reported were associated with blank flow less than the 10% of the smallest reading. Data sampling and system management was performed by a custom-written MATLAB script (The MathWorks Inc., MA), which communicated via serial com ports to the Agilent mass flow meter and pressure controller. This ensured any human variability was minimized and experimental repeatability was maximized.
The nitrogen Q flux through a channel of section A, and length L was modeled as a function of the pressure drop ∆P , as follows: whereP ia the average pressure in the channel, and D is the gas diffusion constant, taken to have Knudsen's form: where h is the channel height, R the universal gas constant, T the temperature, and M the gas molecular weight. This equation is valid only when the gas is confined within a small space where molecule-to-molecule interactions become insignificant compared to molecule-to-surface interaction. This assumption is only true if the ratio between the gas mean free path and the constraining structure are within a certain range. This comparison is quantified by the Knudsen number: where λ is the gas mean free path and L a characteristic length of the system. It is commonly accepted that when K n > 0.05, the system can be described by the Knudsen diffusion equation.
The Knudsen equation was found to reliable predict the gas flow given specified channel geometry and dimensions, as demonstrated by good agreement between the equation and experimental data (see Figure 1j). This test further confirmed the quality of the fabricated membranes.

Supplementary Note 4: Membrane Characterization with Quantum Dot Filtration
CdSe core quantum dots (QDs) (Ocean NanoTech, AR), without substantial surface modification, were solubilized in toluene. QDs were obtained with three different molecular diameters: 2.0, 3.0, and 4.0 nm, which correlated with rated emission peaks at 507, 535, and 584 nm, respectively. The dimensional tolerance was calculated to be ±0.28 nm, which was determined through interpolating the emission peaks for each QD size. Two identical stainless steel reservoirs, each with a 350 µl capacity, were manufactured to tightly fit either side of the nDS membranes (see Supplementary Even with a microfabrication tolerance within 10% of the nanochannel's nominal value, these membranes exhibited QD selectivity with a sub-nanometer resolution of ∼ 0.5 nm. These results strongly evidenced the nanochannels' precise creation.

Supplementary Note 5: Analytes Employed
Histamine (electric charge q = +2e at pH 4) and epinephrine (q = +e at pH 7) were selected as the model cations for this study. Both are hydrophilic, biogenic amine neurotransmitters. Histamine is composed of an imidazole ring and an amino group connected by two methylene groups, which provides an overall estimated molecular diameter of approximately 0.5 nm. Epinephrine is a cathecholamine consisting of a 1,2-dihydroxybenzene connected by hydroxyl and amine groups to a terminating methyl and has an estimated molecular diameter of approximately 0.6 nm. Acetylsalicylic (0 charge at pH 3) acid and L-phenylalanine (0 charge at pH 7) were employed as model neutral solutes. Acetylsalicylic acid (also known as aspirin) is a commonly employed, non-steroidal anti-inflammatory medication. Its structure consists of a benzene ring connected to carboxylic acid and ester functional groups. L-phenylalanine is an essential amino acid with a phenyl, amino, and carboxyl group sequence. Both molecules have an approximate molecular diameter of 0.6 nm.
Cefazolin (q = −e at pH 7) and 3-aminosalicylic acid (q = −e at pH 7) were selected to serve as a model anionic species. Cefazolin is a cephalosporin antibiotic characterized by an elongated conformation and an estimated molecular diameter of approximately 1 nm. 3-aminosalicylic acid is an aminosalicylate employed as an anti-tuberculosis drugs most often with the addition of isoniazid. Its structure consists of a benzene ring connected with amine, hydroxyl, and carboxyl groups and has an approximate molecular diameter of 0.6 nm. Diameters were calculated from estimated molecular volumes [12] and are expected to underestimate the effective hydrodynamic diameters in aqueous solution. Regarding the two neutral molecule employed in the study, aspirin and phenylalanine, it can be shown from Supplementary Figure 2  to drug concentration leveraging the information previously obtained by the standard curves. To achieve the cumulative released drug over time, the latter data was multiplied by the sink volume

Supplementary Note 8: Data Fitting
The diffusivity and permeability of the experimental data were found using the following procedure: • Cumulative release fit with the equation: where n 2 is the total mass released at time t into the sink reservoir, V 1 (200 µ l) and V 2 (4.450 mL) the volumes of the source and sink chambers, respectively, ∆C 0 the initial concentration difference while τ is the equilibration time.
• The equilibration time is related to the membrane permeability P mem [cm 3 s −1 ] by: • Finally, from the permeability, the nanochannel diffusivity D [cm 2 s −1 ] is derived. In the simple case of a single channel, permeability and diffusivity are related by: P mem = whD/L, where w, h and L are the width, height and length, respectively, of the channel.
In our case, the relation between the membrane permeability and the diffusivity in the nanochannels is more complicated, since the permeability contains contributions from the microchannels in the inlet and outlet sections. Taking advantage of the in-series arrangement of the micro-and nanochannels, we can write: where the indices "i", "n" and "o" refer to "inlet", "nanochannels" and "outlet", respectively, N i , N n and N o are the number of micro-or nanochannels in the three respective regions, and D bulk is the bulk diffusivity.
Using this expression, and assuming that the bulk diffusivity is equal to the value D 250 measured in the membranes with the largest nanochannels (h = 250 nm) we extract from the measured permeability the nanochannel diffusivity, D. 3. We compute the electrostatic potential ψ(z) = ψ(z; 0) due to point-charge electrolytes inside the slit-channel; for the sake of simplicity, we consider a slit of height h parallel to the z axis, and infinite in the x and y directions, so that the potential only depends on z; 4. We compute from ψ(z) the surface charge density σ n = 0 k B T dψ(z)/dz| z=h at the position z = h of the silica-diffuse layer interface; 5. We assume that local equilibrium exists between the charges at the channel surface and in the solution, at the silica-diffuse layer interface, and we set σ s = −σ n ; 6. The computed potential depends in general on two parameters: the value of the potential at the silica-diffuse layer interface, ψ d , and the value at the center plane of the channel, ψ n . These two parameters are self-consistently computed by iteratively solving σ s (ψ d , ψ n ) + σ n (ψ d , ψ n ) = 0; 7. We then define the effective diffusivity of the diffusing molecules, as the product of a "repartition coefficient" times the bulk diffusivity. The repartition coefficient accounts for the enhancement, respectively decrease, of the transport rate of positive, respectively negative, solutes interacting with the electrostatic potential created by the surface charge and the NaCl solution within the channel. The repartition coefficient also depends on the steric, hard-sphere interactions of the diffusing molecules with the channel walls.
8. All of the following equations are based on the assumption of two infinitely wide planes and they neglect the effect due to the side walls. Additionally, the reference geometry is as illustrated by the Supplementary Figure 6:

Supplementary Note 10: Silica site-binding model
There are several variants of the silica site-binding model, of varying complexity. In general, these models describe how the silica surface charge varies when exposed to a solution of varying pH or ionic concentration. Assuming local equilibrium, the surface charges have to be neutralized by the charges in solution. If we call σ s the surface charge density, then local equilibrium is guaranteed where σ n is the total charge density in the solution. Following Behrens [1] we can define the diffuse layer potential as a function of the surface charge density starting from the Stern layer potential (ψ 0 ): where σ n is the total charge density in the solution, Θ s is the total concentration of surface sites, 10 pK ≡ K − is the equilibrium dissociation constant of the surface sites and C S is the Stern layer's capacity. As we stated in Section SI 1, the channel walls are not symmetric. However, the model we used is extremely simplified, and taking into account a more realistic asymmetry of the surface charge would only add unnecessary complications, without improving the model outcome. We therefore assumed symmetric and equally negatively charged channel walls throughout.

Supplementary Note 11: Poisson-Boltzmann Equation
The electric double layer (EDL) potential in a channel is usually computed from the Poisson-Boltzmann (PB) equation. Our model assumes a symmetrical solvent (NaCl) and identical silica surfaces on both sides of the nanochannel. Recalling the Gouy-Chapman equation, it is possible to obtain the potential variation for a rectangular nanochannel: Since the surface charge density must balance the charge density in the adjacent solution, it is possible, with the Poisson equation, to relate the surface potential ψ d to the surface charge density σ n with: where ψ d is the Stern diffuse layer electric potential defined in Supplementary Eq. (10), ψ n is the potential at the nanochannel center-line and n 0 is the ionic concentration of the solution. One may note that the effect of pH variation in the bulk can be accounted for by replacing n 0 with n: n = 10 3 N A 10 −pH + 10 pH−14 + n 0 Supplementary Eqs. (10) and (12) allow us to completely specify the electrostatic potential ψ(z) inside the channel.

Supplementary Note 12: Effective Diffusivity
The electrostatic potential ψ(z) dictates the distribution of cations and anions inside the nanochannel through the equilibrium relation : where q = ±Ze, and the Z is the ion's valence.
We now assume that when a charged molecule diffuses though the nanochannel it interacts with an external potential E(z), so that its concentration c(z) also follows a Boltzmann distribution: In particular, we consider E(z) to be the sum of three contributions, E(z) = E HS (z)+E hyd (z)+ qZψ(z). The hard-sphere potential E HS (z) vanishes everywhere inside the channel, and confines the center of mass of the molecules to be at most at a distance r s from the walls at z = 0 and h, where r s is the molecular radius; the term E hyd = k B T ln[K(2r s /h, z/h)], where K −1 (2r s /h, z/h) is the enhanced drag coefficient (see below) results from hydrodynamic interactions.
In conclusion, the molecular concentration inside the channel is where c 0 is an appropriate reference concentration. Note that we are assuming that the diffusing molecules move against a fixed potential and do not modify the charge distribution within the channel. This is clearly not correct, but going beyond this assumption is outside the scope of the paper.
Plecis et al [9] defined a "solute partition coefficient" β, expressed as the ratio between the observed or "effective" permeability P eff , and a "theoretical" permeability P * : Although Plecis et al. interpret the effective permeability as a consequence of a reduced, or effective, channel height, we follow Deen [2] in defining the diffusivity as the quantity that is modified by the confinement inside the channel, so that: Still following Deen [2] , we assume that the partition coefficient is proportional to the average of the solute concentration over the channel height, so that: Supplementary Eq. (16) yields then: The enhanced drag coefficient has been computed for a rectangular channel using the so called "centerline approximation" [2] which consists in letting z = h/2: Finally: In the case of a neutral solute (q = 0) the partition coefficient simplifies to [10,11,2].
A coefficient β(h) > 1 for a given class of molecules implies that their concentration inside the nanochannel is larger than the reference one: this situation is called enrichment. On the contrary, β(h) < 1 for a class of ions implies that those ions are prevented to enter the channel, leading to a lower concentration than in the reference state; this is called exclusion. Supplementary Note 14: Nanochannel pH

Supplementary Note 13: Model Implementation
As described in the previous section, the model iteratively calculates the concentration of ions and, consequently, the concentration of hydronium (H + ) molecules. We were able to observe that the enrichment factor cause an increase H + concentration and, therefore, a reduction of pH.
This phenomenon, in agreement with Plecis [9], is predominant in small channels or at low ionic strength. With this model, we predicted for our smallest nanochannels a decrease in pH to 6.2 from its bulk value 7 (see Supplementary Table 1). The molecules investigated here were chosen in such a way that their physical and chemical properties would not change within such a variation in pH.

Supplementary Note 15: Modeled Diffusivity
This model was compared to experimental diffusivities in Fig. 3 of the main article. Here, we vary the parameters Z and r s in such a way to obtain the best possible fit to the measurements, of the computed diffusivity, with a better quantitative "fit", but still a shape that does match the measurements. Similar considerations apply to neutral molecules, for which this "fit" is even harder to justify.

Supplementary Note 16: Streaming Potential
Measuring the streaming potential allows one to estimate the ζ potential at the silica surfaces of the nanochannels. Applying a mechanical pressure between the nanochannel inlet and outlet forces the molecules to pass through the membrane. The silica surface charges create within the nanochannel an electric field, favoring the flow of counter-ions over that of the co-ions, producing a net ionic current. By doing so, the streaming potential starts to build up between the two floating Ag/AgCl electrodes (see Supplementary Figure 10). At steady state, the potential gradient (U str ) is directly proportional to the ζ potential (ζ) and to the pressure gradient (∆P ) as follows: where r and 0 are respectively the relative and electrical permittivity, η the dynamic viscosity, and K L the specific conductivity of the bulk solution. Employing the Smoluchowski theory as a first approximation, since the thin double layer assumption is not always valid, we can assume ψ (κ −1 ) ≈ ζ where κ −1 is the Debye length. Finally, thanks to the Poisson-Boltzmann equation, the corresponding surface electric potential (ψ 0 ) can be calculated and compared with the one computed from the model.