Non-monotonic variation of flow strength in nanochannels grafted with end-charged polyelectrolyte layers

Abstract

The electrokinetic transport of fluids, also called the electroosmotic flow (EOF), in micro/nanoscale devices occurs in promising applications such as electrokinetic energy conversion (EKEC) systems. Recently, EKEC systems grafted with end-charged polyelectrolyte (PE) layers (PELs) have been reported to exhibit higher efficiencies than those of intrinsic systems. Understanding the interplay between the end-charged PELs and electrical double layers (EDLs) on the EOF is crucial for designing highly efficient EKEC systems. The interplay between the end-charged PELs and EDLs on the strength of the EOF (V₀) is studied by explicitly modeling the EOF through nanochannels grafted with end-charged PELs using atomic simulations. The variation of V₀ is examined for nanochannels grafted with PELs at various separations (d = 3.5–0.4 nm) to cover various conformations of PEs, inlcuding mushroom, semi-dilute brushes, and concentrated brushes. We find that V₀ follows a non-monotonic variation as d decreases and this is correlated with the conformation of the PEs. Specifically, as d decreases, V₀ decreases first in the mushroom regime (d = 3.5–2.0 nm), and then V₀ increases in the concentrated brush regime (d = 0.75–0.4 nm). Navigated by the continuum Navier–Stokes–Brinkman model, the above observations are rationalized by the competition between the driving effect from the spatial shift of ions in EDLs and the drag effect from PELs. The insights obtained in this work are important to guide the design of highly efficient EKEC systems by grafting end-charged PELs onto channel surfaces.

---

1 Introduction

The electrokinetic transport of fluids in micro/nanoscale devices has attracted increasing attention and has been widely used in applications ranging from the sensing and separation of molecules,¹ the gating of liquids,² the gating of ions,³ biomedical care⁴ to energy conversion systems.^5,6 It has been acknowledged that the electroosmotic flow (EOF) can have a significant effect on the transport mechanism of molecules across nanopores.⁷ Electroosmosis can either compete or cooperate with electrophoresis in single molecular trapping in nanopores.⁸ Particular advances have been achieved in energy conversion systems based on electromechanical technologies to harvest energy from renewable sources,⁹ human motion^6,10 or environmental waste heat.¹¹ Electrokinetic energy conversion (EKEC) systems,⁵ which convert energy from micro/nanofluidic flow, have attracted increasing attention owing to their low maintenance and capability to provide a power source at the microscale. The early development of EKEC was pioneered by works from Kwok¹² and Daiguji¹³ with an EKEC energy conversion efficiency of around 1%. With recent advances in nano-materials,¹⁴ nano-manufacturing² and operation parameter optimization (e.g. adjusting temperature¹¹), the EKEC systems have achieved significant progress. In a recent study, an efficiency of up to 50% was achieved for a ballistic electrostatic generator.¹⁵

The key ingredients for the working of EKEC systems include the electrical double layers (EDLs).¹⁶ As a charged or an ionizable substrate is immersed into the electrolyte solution, the charged surfaces attract counter-ions to balance the charge, resulting in the formation of EDLs at the interface of the electrolyte and substrate.^17,18 With an external electric field tangential to the substrate, the ions in the EDLs move and transport the momentum to solvents, which generates the electrokinetic transport of fluids, resulting in electroosmotic flow (EOF).¹⁷ With the absence of the external electric field, the flow of electrolytes in nanochannels can be driven by external pressure differences. The transport of ions results in an electrical current, named the streaming current.¹⁹ The streaming current offers a simple and effective method to convert mechanical energy to electrical power.^12,13,15,20 The electrokinetic transport of fluids plays a crucial role in EKEC devices, which justifies the need for fundamental research on electrokinetic transport.

As the driving force for the electrokinetic transport of fluids occurs at the interface between the substrate and electrolyte solution, the physical and chemical properties of the interface are critical for electrokinetic transport. Grafting polyelectrolyte (PE) brushes to the surfaces is a versatile method to modify the physical and chemical properties of interfaces,^21–23 enabling applications such as ionic gates,²⁴ single (bio)nanoparticle sensing,²⁵ regulating ion transport,²⁶ nanofluidic diodes,²⁷ and current rectification.^28–30 Moreover, the grafting properties of PE brushes may significantly affect the electrokinetic transport in the nanochannel. For example, the Donnan potential^31,32 of the nanochannel significantly affects the electrokinetic transport in a biomimetic PE-modified nanochannel, which is modulated by the imposed gate voltages and the solution properties. The ion selectivity of the biomimetic nanopore and the preferential direction of the ionic current can be regulated by the grafting density of the PE brushes.^23,30 The current rectification and ion concentration polarization effects are strongly affected by the grafting position of the PE brushes at the nanopore (on the inner or outer surfaces of the nanopore).²⁸ Although the EOF in microchannels is significantly suppressed by neutral polymers,^33,34 grafting polyelectrolyte layers (PELs) onto channel surfaces can be an effective way to enhance the charge density of channel surfaces and enhance the EOF velocity.³⁵ Finally, channel surfaces grafted with charged PELs improve the performance of EKEC devices. Das and co-workers^36,37 have reported that the streaming potential in soft channels (channels with PELs) is often larger than that in rigid channels and the efficiency of energy conversion of soft channels is several times larger than that in rigid nanochannels for certain parameters of grafting polymers. Jian et al.¹⁶ reported that, with the combined effects of wall softness and the viscoelastic rheology of fluids, grafting PELs to the channel surfaces improves the efficiency of EKECs under certain optimized parameters. Recently, Das and co-workers^38,39 have reported that grafting end-charged PELs to the channel surfaces greatly improves the strength of the EOF. Different from previous work on the enhancement of flow by increasing the charge density of surfaces,³⁵ the enhancement of flow originates from the spatial shift of ions in the EDLs by the end-charged PELs. In a series of works, Das and co-workers^40,41 have shown that nanochannels grafted with end-charged PELs or poly-zwitterionic layers with certain properties (e.g. grafting density and length of PE) significantly improve the energy conversion efficiency of an EKEC, compared with brush-free nanochannels.

A reliable prediction of electrokinetic transport is required to estimate the optimal operation parameters of EKEC devices, which is crucial for the design of EKEC devices. Although the electrokinetic transport theory has been extensively studied,^42,43 obscure points for electrokinetic transport over bare surfaces, such as flow reversal and the role of the Stern layer, still need to be clarified.⁴⁴ Also, it is challenging to accurately predict the transport of fluids through channels grafted with polymer layers. The classic method is mostly based on the Navier–Stokes (NS) equation for flow, the Poisson–Boltzmann (PB) equation for the distribution of ions, and the Darcy equation for the drag experienced by flow through PELs. However, the method does not accurately model electrokinetic transport through PELs. Firstly, the conformation and hydrodynamic properties of PELs require advanced models, such as the molecular dynamic (MD) simulations and the lattice Boltzmann simulation.⁴⁵ Netz et al.⁴⁶ investigated the effect of electric field on the condensation of PEs and reported the scaling of critical field on the nonequilibrium unfolding of PEs. Secondly, the charged PE beads also affect the distribution of ions in the EDLs. Advanced theories, such as classical density functional theory (cDFT),⁴⁷ are required to properly describe the delicate interplay between charged PE beads and ions. To understand the transport of fluids in PELs, MD simulations have been performed to obtain detailed descriptions of flow transport, distribution of ions, and conformation of polymers at molecular scales. Several groups have applied MD simulations to understand the EOF through polymers. Hickey et al.⁴⁸ studied the EOF through charged polymers and observed that the direction of the EOF reverses well before the net charge of the interface (the wall and PELs) changes sign. Cao et al.⁴⁹ studied the EOF through nanochannels with polymer patterning surfaces by MD simulations and observed that the polymer patterning induced anisotropy of the EOF when the direction of the electric field was changed. In a recent work, Das and co-workers⁵⁰ studied the EOF through nanochannels functionalized with PE brushes by all-atom MD simulations and observed that the direction of the EOF changed by changing the electric field strength.

The above works shed light on the electrokinetic transport of flow in channels grafted with PELs. However, the mechanism for controlling the strength of the EOF is not clear and several key questions remain to be answered. Firstly, in the classic method, the hydrodynamic drag of PELs is accounted for by the friction coefficient governed by Darcy’s equation, and PELs are modeled as resistance centers to flows.³⁸ However, PEs are a string of beads rather than isolated beads and their shielding effect strongly affects their hydrodynamic properties.^45,51 How does a realistic model of PEs affect the variation of the EOF by end-charged PELs? Secondly, as the grafting density of the PEs varies, the conformation of the PEs varies from mushroom to brush-like. How does the conformation of the PEs affect the variation of the EOF by end-charged PELs? Resolving these problems is important to accurately predict the performance of EKEC devices with end-charged PELs. Herein, we study the flow transport through end-charged PELs by MD simulations. We observe that V₀ follows a non-monotonic variation as the separation between PEs (d) decreases. As d decreases from 3.5 nm to 0.4 nm, the conformation of PEs changes from mushroom, to semi-dilute brush, and to concentrated brush. Furthermore, the variation of V₀ strongly correlates with the conformation of the PEs. Specifically, as d decreases, V₀ decreases in the mushroom conformation and then V₀ increases in the concentrated brush conformation.

The rest of the manuscript is organized as follows. The methodology is introduced in Section 2. We propose a technique of velocity decomposition to quantify the competition between the driving effect from ions and the drag effect from PELs. In Section 3, we study flow transport through PELs with various separations and elucidate the mechanism and the structural origin of the variation of flow strength for PELs. We conclude in Section 4.

2 Methods

The EOF in a nanochannel grafted with PELs was simulated by MD simulations. The velocity profiles of the EOF, as well as density profiles of the solvent, ions, and PE beads across the nanochannels, were obtained from the MD simulations. The governing factors on the flow strength were analyzed by the Navier–Stokes–Brinkman (NSB) equation. The NSB equation was applied to compare the competition between the driving effect from net concentration of ions (c_ni) and the drag effect from the PEs.

2.1 MD simulations

The MD system consisted of an electrolyte confined between two parallel walls grafted with two types of polymers (end-charged PELs and neutral polymers). The snapshots of the MD systems are shown in Fig. 1(a) and (b). The ionic strength of the electrolyte (I_bulk) was 3.4 × 10⁻² mol L⁻¹, corresponding to that in experiments (0.56 M).⁵² As the lateral sizes of the systems varied for polymers with various d, the number of ions in the system was different. I_bulk for systems with various d was set to be within a 5% difference. d varied from 0.4 nm to 3.5 nm. For d = 2.0 nm, the grafting density of the polymers (σ) was 0.164 nm⁻², in agreement with that used in experiments.⁵² The walls were constructed by three layers of atoms placed on the FCC lattice with an atom density of 33.3 nm⁻³. For the neutral polymer system, the wall atoms in contact with the electrolyte were charged, with a surface charge density (σ_s) of 3.28 × 10⁻² C m⁻², as used in experiments,⁵² and the atoms of the other layers were left uncharged. For the end-charged PE system, the non-grafted end-beads of PELs uniformly carried charges. The number of charges carried by the end-charged PE beads was the same as that in the neutral polymer system. The charge density of the PELs (σ_s,pel) was equivalent to the surface charge density (σ_s). Due to similar constraints from the sizes of the systems, σ_s,pel for systems with various d was controlled to be within a 5% difference.

Fig. 1 Snapshots from MD simulations and a sketch of the EOF velocity decomposition method. (a and b) A homogeneous electric field E_x was applied in the x direction to drive the EOF through the end-charged PELs grafted to the channel walls. The cations and anions are depicted by red and blue colors, respectively. The solvent is depicted as grey media. The charged wall atoms and end-charged PELs are depicted by the yellow color. The PELs with N = 24 mer were grafted with d = 3.5 nm. (c) A sketch of the decomposition of the EOF velocity V(z) into its components v(z,z′) by the velocity function f(z,z′) for a normalized pulse of c_ni at z′ (i.e., c_ni(z = z′) = 1 and c_ni(z ≠ z′) = 0). The procedure for the decomposition of the EOF velocity is given in the main text. (d) The net concentration of ions c_ni and the velocity function of the flow strength f₀ for PELs with d = 3.5 nm. (e) The components of the flow strength v₀ generated by c_ni. The snapshots of the MD systems were generated by the VMD package.⁶⁹

Depending on the separation between polymers, the lateral dimensions (x and y directions) of the system varied from 5 nm to 7 nm. The channel widths (w) were 27 nm or 37 nm for systems with d > 1.5 nm or d < 1.5 nm, respectively. Such channel widths were wide enough for non-overlapping EDLs in the center of the channel. The strength of the EOF (V₀) was measured from the EOF velocity at the center of the channel. V₀ was independent of w for channels with non-overlapping EDLs.⁵³ w was the distance between the innermost layers of the two walls and z = 0 was chosen at the center of the innermost layer of the bottom wall. A vacuum space with a width of two times the channel width w (∼80 nm) was added in the z direction. The system parameters are summarized in Table S1.†

To focus on the hydrodynamic properties, we used the Weeks–Chandler–Andersen (WCA) potential⁵⁴ for interactions among the solvent, PE beads, and the ions. We chose such a primitive model instead of an atomistic solvent model (e.g. SPC/E model) for the following reasons. First, it is well known that such a primitive solvent model can correctly capture the essential features of the EOF⁵⁴ (e.g. ion distribution across the channel and flow strength of the EOF). Secondly, such a model neglects the chemical details of the solvent molecules and allows us to focus on the hydrodynamic interactions between the solvents and PE brushes. Finally, because of the high thermal noise in the EOF simulations by MD, a large external electric field (above 0.1 V nm⁻¹ (ref. 54–57)) is applied to enhance the statistical accuracy. Such a strong electric field may result in the orientation of solvent molecules and may affect the interfacial properties of fluids, such as the viscosity. We have, therefore, chosen to use such a primitive model by modeling the solvent as non-polar WCA spheres, assuming a background dielectric constant, in order to avoid these complications. The dielectric constant was taken as ε_s = 78 to account for the dielectric properties of the solvent. We note that the permittivity of the solvent in a soft interface may be reduced due to confinement and collective polarization effects.^58,59 However, only the out-of-plane dielectric constant of solvent is reduced, while the in-plane dielectric constant does not change much. We verified by separate simulations that the reduced dielectric constant does not strongly affect the net distribution of ions across the channel. Hence, for simplicity, here we took a fixed dielectric constant across the channel. The force field parameters for the solvent and ions were the same as those in ref. 60. Polymers were modeled by the united-atom model regarding the polyethylene glycol (PEG) molecule. The weighted atoms (carbon and oxygen) in each monomer (mer) of PEG (–CH₂CH₂O–) were represented by WCA spheres. The topology files for the polymers were obtained from the PRODRG2 server.⁶¹ Force field parameters for polymers were taken from the OPLS parameters for hydrocarbons.⁶² Dihedral parameters were tuned to obtain a flexibly deformed polymer. The detailed force field parameters were described in ref. 60.

MD simulations were performed using the simulation package GROMACS 4.5.1⁶². The cut-off radius of the Lennard-Jones (LJ) potential was 2^1/6 σ (σ = 0.3 nm) to model the solution conditions of polymers in a good solvent.⁶⁴ Electrostatic interactions were computed using the particle mesh Ewald (PME) method in 2D (in x and y dimensions)⁶³ with a real-space cut-off radius of 1.3 nm. To generate an EOF, an electric field was applied in the x direction with a strength of 0.08 V nm⁻¹. This field strength was used to enhance the statistical accuracy and lies within the typical range of external electric fields (above 0.1 V nm⁻¹) used in other MD simulations.^49,54–57 In our previous work,⁵¹ we verified that the flow strength of the EOF increases linearly with the external electric field in a range up to 0.16 V nm⁻¹. Therefore, the flow strength in this work holds the linearity with the strength of the external electric field. An NPT simulation of a box of electrolyte with a pressure of one bar was performed to obtain the equilibrium density of the electrolyte in bulk. In the planar system, the density of the electrolyte at the center of the channel is tuned to be the same with that in bulk to control the pressure of the system. The initial configurations of the system were built utilizing Packmol.⁶⁵ The simulations of the planar system were performed in the NVT ensemble for 10 ns to reach a steady state, which were followed by a production run of 100 ns. The time step was 4 fs. The temperature was kept at 300 K by the V-rescaling thermostat.⁶³ Some simulations were also performed with the Nose–Hoover thermostat applied in the orthogonal degrees of the flow direction and quite similar velocity profiles were obtained. In each case, we performed three independent simulation runs with different random seeding. Then, we averaged the results from the independent runs to reduce the fluctuation and obtained error bars for the reported data.

2.2 Navier–Stokes–Brinkman model

The NSB model was adapted from the model in ref. 66.where μ(z) is the fluid viscosity, u_eo is the EOF velocity, c_beads(z) is the number density of the PE beads, a_bead is the effective Stokes radius of the polymer beads, and ϕ_s(z) is the volume of the polymer beads. The function K(ϕ_s(z)) accounts for the correlations between homogeneously distributed spherical particles,⁶⁷ F is the Faraday constant, c_i(z) is the ionic concentration of species i, M is the number of ionic species (M = 2), and E_ext is the applied electric field. On the left-hand side of eqn (1), the first term denotes the viscous force, the second term is the hydrodynamic drag exerted by the polymers on the fluid, and the last term is the driving force due to the external electric field. The driving force of the EOF is the net electrostatic force experienced by the ions (cations and anions) and its magnitude is in proportion to the net concentration of ions c_ni (c_ni = c_cat − c_ani).

The concentrations of ions and PE beads (c_cat(z), c_ani(z), and c_beads(z)) were measured directly in the MD simulations and used as inputs in eqn (1). The fluid viscosity (μ(z)) was obtained from a separate MD simulation using the method described in ref. 60. The model proposed by Batchelor and Green⁶⁸ was used to take into account the variation of the viscosity with the position z across the channel. The no-slip boundary condition was applied at the wall and the mirror symmetry of the EOF velocity profile was imposed at the center of the channel. The hydrodynamic radius of the polymer beads (a_bead), was treated as an adjustable parameter to match the EOF velocity profiles from the continuum model to those obtained by MD simulations.

2.3 EOF velocity decomposition

In light of Green’s function of Poisson’s equation, a velocity function technique was proposed to decompose the EOF velocity into its components generated by c_ni. Due to the linearity of the NSB model (eqn (1)), the velocity profile (V(z)) generated by c_ni across the channel could be decomposed into the components of velocity (v(z,z′)) generated by c_ni at z′. The integral of v(z,z′) over the channel (z′ ∈ (0, H/2)) results in V(z). A schematic sketch of the velocity function technique is shown in Fig. 1(c). The velocity function (f(z,z′)) was the velocity generated by a normalized pulse of c_ni at z′ (i.e., c_ni(z′) = 1 and c_ni(z ≠ z′) = 0). Due to the linearity of the NSB model, v(z,z′) was obtained by scaling f(z,z′) with c_ni(z′), i.e., v(z,z′) = c_ni(z′)f(z,z′).

As for V₀, its velocity function (f₀(z′)) was computed to obtain its components (v₀(z′)). Physically, f₀(z′) quantified the capacity of v₀(z′) generated by a normalized pulse of c_ni at z′. Fig. 1(d) showed that f₀(z′) linearly increased with the position of ions z′ beyond a threshold location, which confirmed the dependency of the capacity of the generated v₀(z′) on the ions position. With calculated f₀(z′), v₀(z′) can be obtained by scaling f₀(z′) with c_ni(z′) (Fig. 1(e)). V₀ is integrated by v₀(z′) (z′ ∈ [0, H/2]). The mathematical formulas of the flow strength decomposition are as follows.

v₀(z′) = f₀(z′)c_ni(z′) (3)

For the same c_ni, we get the same value of V₀ via integration of v₀(z′) and direct solution from the NSB model (shown in Fig. S2(a)†). Besides, the same EOF velocity profile can be derived by these two methods. A detailed derivation of the EOF velocity decomposition is described in Section S2 of the ESI.† f(z,z′) at various z′ are shown in Fig. S2(b)† and V assembled by f(z,z′) is shown in Fig. S2(c).†

3 Results and discussion

Our primary interest is to clarify the mechanism governing the strength of the EOF in the channels grafted with end-charged PELs. As shown in Fig. 1(a) and (b), the distribution of end-charged PE beads results in the spatial shift of ions. The position of end-charged PE beads is related to the conformation of PEs, which is determined by the relative magnitude between the size of the PEs (e.g. the gyration radius, R_g) and d. In the following subsections, we examine the EOF in channels grafted with end-charged PELs with various d to elucidate the factors governing the strength of the EOF.

3.1 The effect of grafting density on the flow strength

Results for the EOF in a channel grafted with end-charged PELs with various d (d = 0.4–3.5 nm) are shown in Fig. 2. As d decreases from 3.5 to 0.4 nm, the conformation of the PEs changes from mushroom, to semi-dilute brush, and to concentrated brush (characterized by a criteria from neutral polymer brushes,²¹ see Section S1 of the ESI† for details). Concentrated PE brushes refer to the self-assembled monolayer (SAM) layers in experiments.⁷⁰

Fig. 2 The EOFs through neutral polymers and end-charged PELs in nanochannels with d = 0.5, 0.75, 2.0, and 3.5 nm from MD simulations. (a) The variation of flow strength through end-charged PELs (V_0,chg) and neutral polymers (V_0,neu) with various d. The blue, green, and pink shadow areas in the plot denote the PEs at mushroom, semi-dilute brushes, and concentrated brushes regimes, respectively. V₀ through the channel with no brushes is shown as a dash line. (b) The concentration of cation (c_cat) and anion (c_ani) across the channel in end-charged PE systems, with the cation shown in solid line and anion shown in dash line. (c and d) The concentration of solvent c_sol (c) and PE beads c_beads (d) across the channel in end-charged PE systems.

The variations of flow strength through end-charged PELs (V_0,chg) and neutral polymers (V_0,neu) are shown in Fig. 2(a). We observe that both V_0,chg and V_0,neu follow a non-monotonic variation. In the mushroom regime (2.0 < d < 3.5 nm), both V_0,chg and V_0,neu decrease as d decreases. In the semi-dilute brush regime (0.75 < d < 2.0 nm), V_0,chg remains almost constant while V_0,neu decreases as d decreases. In the concentrated brush regime (0.4 < d < 0.75 nm), both V_0,chg and V_0,neu increase. The variation of V₀ as a function of d originates from the competition between the spatial shift of c_ni in EDLs and the drag from PELs and will be discussed in the subsequent section. Generally, V_0,chg is larger than V_0,neu over the whole range of d. As shown in Fig. S3,† the ions are partly attracted to the charged walls in the neutral polymer systems, while the ions are attracted to the charged PE beads in the end-charged PE system. Hence, the driving forces in the neutral polymer system partially decreases, resulting in a weaker V₀. In addition, the increase of V_0,neu is larger than that of V_0,chg for polymers with the concentrated brush conformation and V_0,neu approaches V_0,chg at d = 0.4 nm. As shown in Fig. S3(a)–(c),† the peak of the counter-ions reduces for neutral polymers with the concentrated brush conformation, due to the steric effect from the polymers. For polymers with d = 0.4 nm, the distribution of ions between the neutral polymer system and end-charged polymer system is comparable because ions are expelled from the polymers at such grafting density.

To clarify the variation of V_0,chg, we examine the concentration of cations (c_cat) and anions (c_ani) across the channel in end-charged PE systems with various d, as shown in Fig. 2(b). Several separations of PELs (d = 0.5, 0.75, 2.0, and 3.5 nm) are chosen to represent the PEs in the mushroom, semi-dilute brush, and concentrated brush regimes. Compared with those in the neutral polymer system (Fig. S3†), c_cat and c_ani in the end-charged PE system are distributed in a wider region due to the attraction from the end-charged PE beads. As d decreases, cations move away from the walls with an enhanced peak. For PELs with the concentrated brush conformation (d < 0.75 nm), cations form a single peak at the interface between the PELs and the fluid. The distributions of c_cat and c_ani originate from a combination of the electrostatic interaction from the charged PE beads and the steric effect from the PELs. Ions are attracted to the charged PE beads due to electrostatic attraction. Additionally, ions move further away from the wall due to the steric interaction with the PELs.

The distributions of solvents (c_sol) and PE beads (c_beads) across the channel in end-charged PE systems with various d are shown in Fig. 2(c) and (d). c_sol and c_beads are corrolated with the conformation of the PELs. For PELs in the mushroom regime (d > 2.0 nm), c_sol is barely affected by c_beads. For PELs in the concentrated brush regime (d < 0.75 nm), as c_sol in region of PELs (0–5 nm) reduces from ∼30 nm⁻³ to ∼15 nm⁻³, c_beads in the same region of PELs increases from ∼10 nm⁻³ to ∼30 nm⁻³. As d further decreases to 0.4 nm, c_sol in the region of PELs reduces to ∼6 nm⁻³ and c_beads in the region of PELs increases to ∼50 nm⁻³. The repulsion of solvents from the PELs is ascribed to the steric effect from the PE beads for the PELs with concentrated brush conformation.

3.2 The mechanism for the variation of flow strength

To determine the mechanism for the variation of flow strength in end-charged PE systems with different d, the continuum NSB model is applied to quantify the driving and drag effect on the variation of V₀ (Fig. 3). Fig. 3(a) shows the velocity profiles obtained from the MD simulations and those predicted by the NSB model, which shows that the predicted velocity profiles from the NSB model are in good agreement with those from the MD simulations even in the concentrated brush regime, hence the NSB model is used to elucidate the factors governing the flow strength. Other velocity profiles for different d are shown in Fig. S4.†

Fig. 3 Competition between the shift of c_ni and f₀ in end-charged PE systems with various d. (a) Velocity profiles from MD simulations (shown by markers) are compared with those predicted by continuum NSB model (shown by solid lines with the same colors). (b) Velocity function of the EOFs through PELs with various separations d (f_0,d(z)). (c) Distribution of c_ni. (d) The change of the components of flow strength (v₀) from d₁ to d₂ (Δv₀(z)(d₁ → d₂)).

We apply the method of velocity decomposition to quantify the competition between the spatial shift of c_ni and drag from PELs. The procedure of velocity decomposition is described in the Methods section. Essentially, V₀ is decomposed into its components (v₀) perturbed by c_ni at z′. v₀(z) can be quantified through the velocity function (f₀) by scaling f₀(z) with c_ni(z), i.e., v₀(z) = c_ni(z)f₀(z) and the distribution of c_ni for various d (d = 0.5–3.5 nm) are shown in Fig. 3(b) and (c), respectively. For f₀ at the same separation, f₀(z) linearly increases as position z moves away from the walls, which shows that the spatial shift of c_ni away from the walls results in a stronger driving effect. For the velocity function with different d (f_0,d), f_0,d decreases as d decreases due to the increased drag. For the distribution of c_ni in end-charged PE systems with various d (Fig. 3(c)), they show distinct features depending on the conformation of the PELs. Specifically, c_ni remains almost constant for PELs in the mushroom regime (d = 2.0–3.5 nm) and it shifts away from the walls by a larger magnitude for PELs in the brush-like regime (d = 0.5–2.0 nm).

To elucidate the competition between the spatial shift of c_ni and drag from the PELs, the change of v₀ (Δv₀(d₁ → d₂)) as d varies is plotted in Fig. 3(d). Specifically, Δv₀(z)(d₁ → d₂) denotes the change of v₀(z) as d decreases from d₁ to d₂, i.e., Δv₀(d₁ → d₂) = v₀(d₁) − v₀(d₂). The positive part of Δv₀(z)(d₁ → d₂) results in a decrease of V₀ (V_0,d1 > V_0,d2) and the negative part of Δv₀(z)(d₁ → d₂) results in an increase of V₀ (V_0,d1 < V_0,d2). The change of v₀(z) shows unique features depending on the conformation of the PELs. In the mushroom regime, the positive part of Δv₀(z)(3.5 → 2.0) spans from 2 nm to 3 nm and its negative part is almost zero, which shows that the drag from the PELs increases by a larger magnitude than the spatial shift of c_ni, as supported by a constant position of c_ni for d = 3.5–2.0 nm. In the semi-dilute brush regime, the positive part of Δv₀(2.0 → 0.75) is comparable to its negative part, which implies that the increase of drag is comparable to the spatial shift of c_ni, in line with a moderate shift of c_ni for d = 2.0–0.75 nm. In the concentrated brush regime, the positive part of Δv₀(0.75 → 0.5) is smaller than the negative part, which indicates that the drag from the PELs increases by a smaller magnitude than the spatial shift of c_ni, as supported by a significant shift of c_ni with an enhanced peak for d = 0.75–0.5 nm.

3.3 Structural analysis of PELs with various separation

The aforementioned analysis shows that the spatial shift of c_ni correlates with the variation of V₀. Meanwhile, c_ni is affected by the conformation of the PELs due to the electrostatic attraction from the end-charged PE beads. Therefore, we analyze the structure of the PELs to elaborate upon the structural origin of c_ni in Fig. 4.

Fig. 4 Structural properties of the PELs with various d. (a) Distribution of the charges of the PELs (ρ_chg). (b) Probability distribution of the parallel orientation angle (P(θ_x)). (c and d) Net radial number density of ions (Δρ_i(r)) (c) and radial number density of solvents (ρ_sol(r)) (d) around end-charged PE beads.

Fig. 4(a) shows the distribution of the charges of the PELs (ρ_chg) in systems with various d. The distribution of ρ_chg is calculated by the volume charge density of the end-charged PELs. The distribution of ρ_chg shows similar features to c_ni. Therefore, the distribution of ρ_chg is analyzed to elaborate the structural origin of c_ni. As d decreases, the spatial shift of ρ_chg increases as the space available for the polymer beads to pack decreases. Specifically, ρ_chg shifts by a small amount in the mushroom conformation since there is enough space for the PE beads to pack. The spatial shift of ρ_chg increases in the brush-like conformation. ρ_chg shift by a significant amount in the concentrated brushes conformation. This can be ascribed to the dense layer formed by the polymers where ions and solvents get expelled from PELs. Under such confinement the PE beads have insufficient space to pack in the polymers layer.

To elucidate the arrangement of PEs, the distribution of the parallel orientation angle (θ_x) is calculated in Fig. 4(b), where θ_x is the angle formed between the orientation vector of the PEs and the unit vector of the x axis. The orientation vector of the PEs points from the grafted bead at the wall to the end-charged bead. As d decreases, the probability of θ_x (P(θ_x)) shifts from 45° to 90°, corresponding to the unfolding of PEs from the mushroom to the brush-like conformation. For polymers with d = 3.5 nm, P(θ_x) at 45° is higher than that at 135°, due to the inducement of flow along the x axis. The probability of the perpendicular orientation angle (θ_z) (P(θ_z)) changes from 45° to 0° (shown in Fig. S5(a)†), corresponding to the unfolding of the PEs from the mushroom to brush-like conformation.

The solvation of end-charged beads by species is quantified by the net radial number density of ions(Δρ_i(r)), the radial number density of PE beads (ρ_beads(r)) and solvents (ρ_sol(r)) around the end-charged PE beads, as shown in Fig. 4(c) and (d) and S5b.† Δρ_i(r) and ρ_sol(r) show different features which are strongly dependent on the conformation of the PEs. In the mushroom regime, as d decreases, Δρ_i(r) reduces by a large amount (the value of first peak of Δρ_i(r) reduces by a factor of ∼10), while ρ_sol(r) slightly reduces. This is beacuse that as d decreases, the amount of end-charged PE beads increases by a factor of ∼3, while the amount of c_ni remains constant, resulting in a significant reduction of Δρ_i(r) in the mushroom conformation. A sufficient amount of solvents surrounds the end-charged PE beads in the mushroom conformation, resulting in a slight reduction of ρ_sol(r) in such conformation. The variation of Δρ_i(r) and ρ_sol(r) in the mushroom regime shows that the space in the PELs in such regime is large enough that solvents and ions are not expelled from the PELs. In the concentrated brush conformation, as d decreases, Δρ_i(r) decreases by a small amount, while ρ_sol(r) is reduced by a large amount (the value of first peak of ρ_sol(r) reduces by ∼50%). This is because that as d decreases, the amount of charged PE beads increases by a factor of ∼2, while the amount of ions around the end-charged PE beads also increases due to the enriched ions expelled from PELs, resulting in a weak decrease of Δρ_i(r) in the concentrated brush conformation. Solvents are expelled from the PELs and the amount of solvents around the end-charged PE beads reduces, resulting in a large reduction of ρ_sol(r) in the concentrated brush conformation. The variation of Δρ_i(r) and ρ_sol(r) in the concentrated brush conformation shows that the space in the PELs with such conformation is so limited that species are expelled from the PELs. The evolution of the solvation structure of the end-charged beads is consistent with the spatial shift of c_ni and they are all correlated with the conformation of PEs.

3.4 Other factors affecting the variation of flow strength

3.4.1 Flow enhancement by end-charged PELs. In the MD systems, V₀ in channels grafted with PELs is weaker than that with no brushes, which is different from other studies.^38,40,41 This may be due to the short polymer length that can be modeled by the current MD simulations. The polymer lengths (N = 2000 mer) in other studies^38,40,41 are much larger than that in the current system (N = 24 mer).

3.4.2 The effect of PELs with fixed end-charged PE beads on the flow strength. In continuum models,^16,38 the PELs are usually modeled as fixed polymer beads. While in MD simulations, PELs are modeled as freely moving polymers. V₀ through PELs with fixed end-charged PE beads (fixed-PELs) and freely moving end-charged PE beads (freely-PELs) are shown in Fig. S6.† The end-charged PE beads in the fixed-PEL system are at the peak position of ρ_chg in the freely-PEL system. A comparison of c_beads between fixed-PEL system and freely-PEL system is shown in Fig. S7.† Fig. S6† shows that V₀ through fixed-PELs is stronger than that through freely-PELs in the brush-like regime, which originates from a more concentrated distribution of c_ni at the fixed-PELs system, leading to a stronger driving effect. In the mushroom regime, V₀ through fixed-PELs is comparable with that through freely-PELs. This behavior may be caused by the fact that ions are less affected by the position of charged PE beads because there is more space available in PELs with the mushroom conformation.

3.4.3 The effect of the space charge density of PE brushes on the flow strength. In fact, the charged functional groups may not be at the tail of a PE brush. We studied the effect of charged polymer beads distributed at the middle and in the tail of the PE brushes. The results show that V₀ through such PE brushes is generally weaker than that through the end-charged PE brushes (see Section S3 of the ESI† for details).

3.4.4 The effect of the ionic strength of the solution on the flow strength. The ionic strength of the electrolyte solution is another factor affecting the variation of flow strength. We examine the effect of ionic strength on the flow strength, as shown in Fig. S9.† Due to the computational limit, the effect of ionic strength is examined by increasing the ionic strength I by a factor of 3.5, 5, and 7.5. The results show that V₀ decreases as the ionic strength increases for PELs in the mushroom conformation. Such an observation agrees with previous work for cases with short polymer lengths.³⁸ In the current system, as the ionic strength increases, c_ni shifts towards the walls, resulting in a weaker driving effect of c_ni and the decrease of V₀.

3.4.5 Hydrodynamic radius of the PE beads. The hydrodynamic radius of the PE beads (a_bead) is the only fitting parameter used by the NSB model to match the velocity profiles of MD simulations. We examine the a_bead of end-charged PELs with various d in Fig. S10.† a_bead is around 0.1 of its physical size a₀ (i.e., LJ size and a₀ = 0.156 nm) for PELs with the mushroom conformation, in agreement with our previous work.⁵¹ As d decreases, a_bead increases until ∼0.5a₀ for PELs with the concentrated brush conformation. The increase of a_bead may be due to the unfolding of PEs from mushroom to brush-like configuration. From our previous work,⁵¹ a_bead in the brush-like configuration is larger than that in mushroom configuration.

4 Conclusions

We have studied the variation of V₀ in nanochannels grafted with end-charged PELs using MD simulations, analyzed by the continuum NSB model. We observe that V₀ follows a non-monotonic variation as d decreases from 3.5 nm to 0.4 nm. V₀ decreases first as 2.0 < d < 3.5 nm, in agreement with previous studies.³⁸ In particular, V₀ increases as d decreases for PELs with a concentrated brush conformation (d < 0.75 nm), a new behavior of flow strength.

The variation of V₀ through PELs with d results from the competition between the increased drag from PELs and the spatial shift of c_ni away from the walls. A method using the velocity function (f₀) is proposed to decompose V₀ into v₀ (the components of V₀ perturbed by c_ni). Δv₀ (the change of v₀) is examined to quantify the competition between the driving effect from c_ni and the drag effect from PELs. For PELs in the mushroom regime (2.0 < d < 3.5 nm), the drag from PELs dominates over the driving effect from the spatial shift of c_ni, resulting in the decrease of V₀. For PELs in the concentrated brush regime (d < 0.75 nm), the driving effect from the spatial shift of c_ni dominates over the drag from PELs, resulting in the increase of V₀. The structural analysis of the PELs shows that ρ_chg is highly correlated with the distribution of c_ni. The distribution of ρ_chg is rationalized by the conformation of the PELs. In the mushroom regime, there is enough space in the PELs for the PE beads to pack, therefore ρ_chg changes little as d decreases. In the concentrated brush conformation, less space is available in the PELs for the PE beads to pack, resulting in the expulsion of ions and solvents from the PELs, therefore ρ_chg shifts by a large amount as d decreases. The radial density distribution of species also demonstrates that enough spaces are available in the PELs with the mushroom conformation for ions and solvents to pack, leading to reduced ions around the end-charged PE beads. While there is less space available in the PELs with the brush conformation, resulting in the expulsion of ions from the PELs and enhanced ions around the end-charged PE beads. Overall, the variation of V₀ through end-charged PELs originates from the interplay between the structure of the PELs and the distribution of ions. Although the enhancement of the electrokinetic transport over concentrated brush-like PELs has not been validated by experimental measurements, it is well known that high packing densities are crucial for the tribological performance of monolayers.^71–73 Our theoretical work suggests that densely packed monolayers can enhance electrokinetic transport and encourages further experimental studies.

Author contributions

P. W.: conceptualization, funding acquisition, supervision and writing – original draft; T. S.: software, data curation, visualization, investigation and writing – review & editing; X. J.: investigation and writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

P. W. gratefully acknowledges the support from the Natural Science Foundation of Hubei Province of China under grant No. 0216120087. The authors thank Prof. Svyatoslav Kondrat and Dr Ming Chen for reading the manuscript and providing numerous suggestions.

Notes and references

1. Z. Adamczyk, B. Jachimska, T. Jasiński, P. Warszyński and M. Wasilewska, Colloids Surf., A, 2009, 343, 96–103.
2. J. Zhang, K. Zhan, S. Wang and X. Hou, Soft Matter, 2020, 16, 2915–2927.
3. Z. Adamczyk, K. Sadlej, E. Wajnryb, M. L. Ekiel-Jezewska and P. Warszynski, J. Colloid Interface Sci., 2010, 347, 192–201.
4. R. Qiao and P. He, Langmuir, 2007, 23, 5810–5816.
5. S. Das, M. Banik, G. Chen, S. Sinha and R. Mukherjee, Soft Matter, 2015, 11, 8550–8583.
6. K. Liu, T. Ding, J. Li, Q. Chen, G. Xue, P. Yang, M. Xu, Z. L. Wang and J. Zhou, Adv. Energy Mater., 2018, 8, 1702481.
7. G. Huang, K. Willems, M. Soskine, C. Wloka and G. Maglia, Nat. Commun., 2017, 8, 935.
8. M. Chinappi, M. Yamaji, R. Kawano and F. Cecconi, ACS Nano, 2020, 14, 15816–15828.
9. A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell and L. Bocquet, Nature, 2013, 494, 455–458.
10. K. Liu, T. Ding, X. Mo, Q. Chen, P. Yang, J. Li, W. Xie, Y. Zhou and J. Zhou, Nano Energy, 2016, 30, 684–690.
11. R. Long, F. Wu, X. Chen, Z. Liu and W. Liu, Int. J. Heat Mass Transfer, 2021, 168, 120842.
12. J. Yang, F. Lu, L. W. Kostiuk and D. Y. Kwok, J. Micromech. Microeng., 2003, 13, 963–970.
13. H. Daiguji, P. D. Yang, A. J. Szeri and A. Majumdar, Nano Lett., 2004, 4, 2315–2321.
14. R. Qu, X. Zeng, L. Lin, G. Zhang, F. Liu, C. Wang, S. Ma, C. Liu, H. Miao and L. Cao, ACS Nano, 2020, 14, 16654–16662.
15. Y. Xie, D. Bos, L. J. de Vreede, H. L. de Boer, M.-J. van der Meulen, M. Versluis, A. J. Sprenkels, A. van den Berg and J. C. T. Eijkel, Nat. Commun., 2014, 5, 3575.
16. Y. J. Jian, F. Q. Li, Y. B. Liu, L. Chang, Q. S. Liu and L. G. Yang, Colloids Surf., B, 2017, 156, 405–413.
17. J. Lyklema, Fundamentals of Interface and Colloid Science, Academic Press, San Diego, CA, 1995.
18. R. J. Hunter, Zeta Potential in Colloid Science: Principles and Applications, Academic Press, London, 1981.
19. F. H. J. van der Heyden, D. Stein and C. Dekker, Phys. Rev. Lett., 2005, 95, 116104.
20. C. Davidson and X. Xuan, J. Power Sources, 2008, 179, 297–300.
21. W. L. Chen, R. Cordero, H. Tran and C. K. Ober, Macromolecules, 2017, 50, 4089–4113.
22. X. Hou, W. Guo and L. Jiang, Chem. Soc. Rev., 2011, 40, 2385–2401.
23. Z. Zeng, L.-H. Yeh, M. Zhang and S. Qian, Nanoscale, 2015, 7, 17020–17029.
24. Z. Zeng, Y. Ai and S. Qian, Phys. Chem. Chem. Phys., 2014, 16, 2465–2474.
25. M. Ali, S. Nasir, P. Ramirez, J. Cervera, S. Mafe and W. Ensinger, J. Phys. Chem. C, 2013, 117, 18234–18242.
26. L.-H. Yeh, M. Zhang, S. Qian, J.-P. Hsu and S. Tseng, J. Phys. Chem. C, 2012, 116, 8672–8677.
27. P.-H. Peng, H.-C. Ou Yang, P.-C. Tsai and L.-H. Yeh, ACS Appl. Mater. Interfaces, 2020, 12, 17139–17146.
28. J.-P. Hsu, H.-H. Wu, C.-Y. Lin and S. Tseng, Phys. Chem. Chem. Phys., 2017, 19, 5351–5360.
29. C. Y. Lin, J. P. Hsu and L. H. Yeh, Sens. Actuators, B, 2018, 258, 1223–1229.
30. T.-W. Lin, J.-P. Hsu, C.-Y. Lin and S. Tseng, J. Phys. Chem. C, 2019, 123, 12437–12443.
31. L. Benson, L.-H. Yeh, T.-H. Chou and S. Qian, Soft Matter, 2013, 9, 9767–9773.
32. Z. Milne, L.-H. Yeh, T.-H. Chou and S. Qian, J. Phys. Chem. C, 2014, 118, 19806–19813.
33. M. Monteferrante, S. Melchionna, U. M. B. Marconi, M. Cretich, M. Chiari and L. Sola, Microfluid. Nanofluid., 2015, 18, 475–482.
34. M. Monteferrante, L. Sola, M. Cretich, M. Chiari, U. Marini Bettolo Marconi and S. Melchionna, J. Chem. Phys., 2015, 143, 184907.
35. J. Li and D. Q. Li, J. Colloid Interface Sci., 2019, 553, 31–39.
36. S. Chanda, S. Sinha and S. Das, Soft Matter, 2014, 10, 7558–7568.
37. J. Patwary, G. Chen and S. Das, Microfluid. Nanofluid., 2016, 20, 37.
38. G. Chen and S. Das, J. Phys. Chem. B, 2017, 121, 3130–3141.
39. V. S. Sivasankar, S. A. Etha, H. S. Sachar and S. Das, Phys. Rev. E, 2020, 102, 013103.
40. G. Chen, H. S. Sachar and S. Das, Soft Matter, 2018, 14, 5246–5255.
41. G. Chen, J. Patwary, H. S. Sachar and S. Das, Microfluid. Nanofluid., 2018, 22, 112.
42. W. B. Russel, D. A. Saville and W. R. Schowalter, Colloidal Dispersions, Cambridge University Press, Cambridge, 1989.
43. P. H. Wiersema, A. L. Loeb and J. T. G. Overbeek, J. Colloid Interface Sci., 1966, 22, 78–99.
44. M. Rezaei, A. R. Azimian and D. T. Semiromi, Heat Mass Transfer, 2015, 51, 661–670.
45. S. Raafatnia, O. A. Hickey and C. Holm, Macromolecules, 2015, 48, 775–787.
46. R. R. Netz, Phys. Rev. Lett., 2003, 90, 12.
47. X. Hu, X. Kong, D. Lu and J. Wu, J. Chem. Phys., 2018, 148, 084701.
48. O. A. Hickey, C. Holm, J. L. Harden and G. W. Slater, Macromolecules, 2011, 44, 9455–9463.
49. Q. Cao, Soft Matter, 2019, 15, 4132–4145.
50. T. H. Pial, H. S. Sachar, P. R. Desai and S. Das, ACS Nano, 2021, 15, 6507–6516.
51. P. Wu, T. Sun, X. Jiang and S. Kondrat, Polymers, 2019, 11, 1038.
52. J. A. Cohen and V. A. Khorosheva, Colloids Surf., A, 2001, 195, 113–127.
53. R. Qiao and N. R. Aluru, Int. J. Multiscale Comput. Eng., 2004, 2, 173.
54. F. Tessier and G. W. Slater, Macromolecules, 2006, 39, 1250–1260.
55. R. Qiao and N. R. Aluru, J. Chem. Phys., 2003, 118, 4692–4701.
56. R. Qiao and N. R. Aluru, Colloids Surf., A, 2005, 267, 103–109.
57. R. Qiao and N. R. Aluru, Langmuir, 2005, 21, 8972–8977.
58. A. Schlaich, E. W. Knapp and R. R. Netz, Phys. Rev. Lett., 2016, 117, 048001.
59. L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Janardanan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe, G. Gomila, K. S. Novoselov and A. K. Geim, Science, 2018, 360, 1339–1342.
60. P. Wu and R. Qiao, Phys. Fluids, 2011, 23, 072005.
61. A. W. Schüttelkopf and D. M. F. van Aalten, Acta Crystallogr., Sect. D: Biol. Crystallogr., 2004, 60, 1355–1363.
62. S. W. I. Siu, K. Pluhackova and R. A. Boeckmann, J. Chem. Theory Comput., 2012, 8, 1459–1470.
63. E. Lindahl, B. Hess and D. van der Spoel, J. Mol. Model., 2001, 7, 306–317.
64. G. S. Grest, Phys. Rev. Lett., 1996, 76, 4979–4982.
65. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
66. R. J. Hill, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 051406.
67. D. L. Koch and A. S. Sangani, J. Fluid Mech., 1999, 400, 229–263.
68. G. K. Batchelor and J. T. Green, J. Fluid Mech., 1972, 56, 401–427.
69. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
70. D. T. S. Ranathunga, A. Shamir, X. Dai and S. O. Nielsen, Langmuir, 2020, 36, 7383–7391.
71. J. E. Black, A. Z. Summers, C. R. Iacovella, P. T. Cummings and C. McCabe, Nanomaterials, 2019, 9, 639.
72. B. D. Booth, S. G. Vilt, C. McCabe and G. K. Jennings, Langmuir, 2009, 25, 9995–10001.
73. B. Bhushan, T. Kasai, G. Kulik, L. Barbieri and P. Hoffmann, Ultramicroscopy, 2005, 105, 176–188.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1ra06601c

---