Electro-osmotic diode based on colloidal nano-valves between double membranes

The rectification of electro-osmotic flows is important in micro/nano fluidics applications such as micro-pumps and energy conversion devices. Here, we propose a simple electro-osmotic diode in which colloidal particles are contained between two parallel membranes with different pore densities. While the flow in the forward direction just pushes the colloidal particles toward the high-pore-density membrane, the backward flow is blocked by the particles near the low-pore-density membrane, which clog the pores. Nonequilibrium molecular dynamics simulations show a strong nonlinear dependence on the electric field for both the electric current and electro-osmotic flow, indicating diode characteristics. A mathematical model to reproduce the electro-osmotic diode behavior is constructed, introducing an effective pore diameter as a model for pores clogged by the colloidal particles. Good agreement is obtained between the proposed model with estimated parameter values and the results of direct molecular dynamics simulations. The proposed electro-osmotic diode has potential application in downsized microfluidic pumps, e.g., the pump induced under AC electric fields.


INTRODUCTION
Micro-and nano-scale transport systems have attracted interest for application in fields ranging from bioand medical-technologies to new energy conversion and storage [1][2][3][4][5][6][7][8][9]. Various devices exploiting the phenomena particular to the small scale flows have been proposed [10][11][12]. In typical micro-and nano-fluidic devices, an external field, such as an electrical or chemical potential field, applied to an electric double layer formed in the vicinity of a solid-liquid interface, induces electrokinetic phenomena, which are typically accompanied by interfacial liquid flow. The high surface area to volume ratio of micro-and nano-scale systems enhances these phenomena, which have thus been studied for application as the driving force in small-scale transport systems [13,14].
A representative electrokinetic phenmenon is the electro-osmotic flow, which is induced by an external electrical field [15][16][17][18][19]. Since the flow is directly induced by the electric field, i.e., there are no moving parts, it is recognized as a key technique for downsizing of pumps. Applications of porous materials, e.g. a glass porous material, as a source of many pores were extensively studied as a promising setup for downsized pumps [20][21][22]. The use of small size fabrication techniques such as microlithography and chemical etching techniques on a substrate was alternatively studied [23,24]. Recently, all-plastic nano fabrication of electro-osmotic flow membrane was reported [25]. Along with these experiments, theoretical analyses for predicting the flow rate [26][27][28][29] and molecular dynamics (MD) simulations [30][31][32][33] have also been conducted for nano sized systems.
One of the challenges for the practical application of electro-osmotic flow is that the long-term application of a DC voltage as an external field electrolyzes water. To prevent the electrolysis of solvent water, many attempts have been made by using, e.g., traveling-wave poten-tials [34], three-dimensional stepped electrode arrays [35], ratcheted electrodes [36], and an AC voltage whose period is shorter than a characteristic time scale of electrolysis [37,38]. In these methods, rectifying the forward and backward electro-osmotic flows is commonly a key ingredient, to induce a net one-way flow.
In the present study, we propose an electro-osmotic diode made of two electro-osmotic flow membranes, that quasi-mechanically suppress backward flow when an electric field is applied in the backward direction. More specifically, the system consists of two membranes with different pore densities and colloidal particles in between (see Fig. 1(a)). When an electric field is applied in the forward direction, the colloidal particles approach the high-pore-density membrane, which has little impact on the generated electro-osmotic flow. The other membrane has much fewer pores, so when an electric field is applied in the opposite direction, the colloidal particles plug the pores, blocking the backward electro-osmotic flow. This simple diode system highlights potential application in creating a one-way flow using the electro-osmotic flow under an AC electric field as mentioned above, if the typical period is sufficiently longer than the relaxation time of the colloid motion.
We conduct MD simulations of the proposed system to demonstrate that both the electric current and electroosmotic flow rate exhibit strong nonlinearity with respect to the electric field, to show a performance as an electro-osmotic diode. In order to confirm the observation, next we construct a mathematical model of the electro-osmotic diode. Here we introduce a new parameter, what we call an effective pore diameter, which changes depending on the existence probability of colloidal particles near pores. The theoretical model for electro-osmotic flow through a cylindrical pore of a finite length is extended using the effective pore diameter. The existence probability density functions of colloidal particles near pores are expressed in terms of the solutions of the Fokker-Planck equations for the motion of colloidal particles. The parameters in the model equation are estimated using independent MD simulations, which are dedicated to measure the properties of colloids. The proposed model with these estimated parameters is shown to reproduce the simulation results for the entire system.

Electro-osmotic diode
Let us consider a system of two membranes M 1 and M 2 , placed in parallel in the direction perpendicular to the z axis, and colloidal particles between the membranes. The particle size is slightly larger than the pore diameter, D, in the membranes. As illustrated schematically in Fig. 1(a), the number of pores for membrane M 1 is comparable to that of the colloidal particles, whereas that for membrane M 2 is higher. The colloidal particles are positively charged, with a charge q per particle, and the walls of the pores are negatively charged. If an electric field E z in the z direction is applied from M 1 to M 2 , the colloidal particles are induced to flow by electrophoresis toward membrane M 2 . An electric field applied in the opposite direction would drive the colloidal particles toward M 1 . In the latter case, we expect that the current and electro-osmotic flow induced by E z (< 0) will be blocked by particles near pores.
The systems considered here for the simulations and theories are on the scale of nm ∼ tens of nm. Corresponding experimental setups could be constructed using materials such as polycarbonate [25], PET [37,79], and car- and (b) the model system considered in the simulations. In the model system, membrane M1 has a physical pore, and membrane M2 is modeled as an energy barrier that affects only the colloid particle.
bon [80,81] for the membranes. The track-etching technique [82,83] is able to control the size and the number of pores, as done for creating nanofilteration membranes for molecular sieving. The electro-osmotic flow is then created using the aqueous electrolyte solution as NaCl solution. Nano-particles such as gold nanoparticles [84], silica particles [85], and hydrous zirconia particles [86] are possible candidates for the colloidal particles.
In the present study, to focus on the interaction between the pore in membrane M 1 and the colloidal particles, the model system shown in Fig. 1(b) is considered as an abstracted system for analysis, where membrane M 1 with thickness L has a single pore with a diameter D, and a single colloidal particle exists in the solution. Membrane M 2 is modeled as a completely semi-permeable membrane with a virtual energy barrier that affects only colloidal particles, i.e., the colloidal particles are completely rejected by M 2 whereas the solvent freely passes through. Assuming periodic boundary conditions in the x and y directions, we investigate the behavior of the current density and electro-osmotic flow of this system using MD simulations.

Molecular dynamics simulations
In this study, MD simulations are performed with the membranes, colloidal particle, and solvent modeled using particles. The Lennard-Jones (LJ) potential is used for inter-particle interactions. The potential energy for the pairwise interaction is E = 4 ((σ/r * ) 12 − (σ/r * ) 6 )) for r * < r c , and zero otherwise. Here, r * is the distance between the center of the interacting particles, and and σ are parameters corresponding to the potential depth and diameter of the particles, respectively; r c is the cutoff parameter. In the present study, common values of LJ parameters, namely 0 and σ 0 , are assigned to all the LJ particles, and the cutoff distance is r c = 2.5σ 0 .
A snapshot of the MD simulation, illustrating the computational system, is shown in Fig. 1(b). We use a sevenlayer hexagonal close-packed structure stacked in the z direction as membrane M 1 , the middle layer of which is placed at z = L z /4, with L z being the system size in the z direction. A pore at the center of the membrane is created by removing particles less than 1.2σ 0 away from the center line. The uniformly distributing 12 particles on the pore wall are monovalent anions, where the charge of each anion is −q 0 . Membrane M 2 , which is a semipermeable membrane, is modeled by a virtual energy barrier that affects only the colloidal particle. Specifically, we assume the following Gaussian potential barrier: where the height of barrier U 0 is set to 30 0 , which is sufficiently large to prevent the colloidal particle from passing through M 1 . The center of the barrier is at z 0 = 3L z /4 and the parameter a 0 , which determines the potential width, is 10σ −2 0 . The colloidal particle is composed of 12 cations located at the vertices of an icosahedron with a diagonal length of 2 σ 0 . The net charge of the colloidal particle is q, i.e., the charge of each ion is q/12. The electrolyte solution consists of 2410 particles, including 20 monovalent anions, and 32 − q monovalent cations, where the number of cations is determined such that the total charge in the system (including the surface charge on the pore wall) is zero.
We use the open-source package LAMMPS [87] for the MD simulations. The velocity Verlet method is employed for the time integration of the Newton equation for each particle. To deal with the long-range Coulomb interaction, we employ the particle-particle-particle-mesh (PPPM) method [88]. A periodic boundary condition is applied in all directions and the NVT ensemble is used to maintain the temperature of the system at T 0 . The time step is dt = 0.002 τ 0 . The typical size of the system is L x = L y = 10 σ 0 in the x and y directions and L z is L z ∼ 40 σ 0 . The precise value of L z is determined such that the pressure in the z direction is 1.0 ± 0.1 0 /σ 3 0 .
The simulation results for various values of colloidal particle charge q are shown in Fig. 2. The current I z and the electro-osmotic flow rate Q z are plotted as functions of the electrical field E z in the range −0.7E 0 ≤ E z ≤ 0.7E 0 . To obtain these results, the system is first equilibrated with no electric field for 10 5 time steps before the production run, with the pressure kept at 0 /σ 3 0 and the temperature kept at T 0 . Then the electric field is turned on and the production run is carried out for more than 10 7 steps. The current and flow rate at a time instance are calculated as I z = (S/V ) ions q i v zi and Q z = (1/N ) solution v zi , where q i and v zi are the charge and velocity in the z direction of a particle, S and V are respectively the area in the x-y plane and the volume of the measured region, and N is the number of the measured particles. These instant values are averaged over the last 8×10 6 time steps of the production run. We perform three production runs for each case using different initial configurations. Each point in Fig. 2 is the mean of three values and the lines are the results of the model (see Sec. for details).
Both the current I z and flow rate Q z for q > 0 clearly exhibit a nonlinear dependence on the electric field E z . Generally, the flow for E z < 0, the backward flow, is suppressed, and completely blocked for q = 6q 0 . For E z < 0, the colloidal particle is dragged close to the pore of M 1 and serves as a valve that blocks the current and flow through the pore. A large colloidal particle charge enhances this effect. To demonstrate the behavior of the colloidal particle, Fig. 3 shows the probability distribution function for the colloidal particle along the z axis for q = 6q 0 and E z = ±0.7E 0 . It is confirmed that the colloidal particle is mostly near membrane M 1 under the backward electric field (E z = −0.7E 0 ), whereas it is sufficiently far from M 1 under the forward electric field (E z = 0.7E 0 ).
The MD results suggest that the nonlinear behavior in the current and flow stems from the effective pore size being decreased by the colloidal particle. We thus construct a theoretical model that captures the nonlinear responses, introducing an effective pore diameter that depends on the existence probability of the colloidal particle in the vicinity of the pore. Using the effective pore diameter, we extend the theoretical model proposed by Sherwood et al. [27], which gives the current and electroosmotic flow induced in a cylindrical pore of finite length. The model equations are given in the next section, followed by a comparison between the model and the results of the MD simulations in Sec. .

MODELING OF ELECTRO-OSMOTIC DIODE
In this section, we develop a mathematical model for the electro-osmotic diode described in the previous section to reproduce the MD results of Fig. 2 outline the model equations for the current and electroosmotic flow in a cylindrical pore of finite length proposed by Sherwood et al. [27], which we employ to express flows without colloidal particles. We then develop a model equation for the effective pore diameter, which depends on the existence probabilities of colloidal particles near the pore. The Fokker-Planck equations governing the existence probabilities are summarized and the analytical solutions are derived.

Electro-osmotic flow through nano-pore
We begin the model construction by considering the current and electro-osmotic flow through a cylindrical pore of finite length L and diameter D as shown in Fig. 4. Here, we give a brief derivation of the model equations, the details of which are found in Refs. [27,28].  In the linear response regime, the current I z and electro-osmotic flow rate Q z are written as follows: where G is the conductance, H is the electro-osmotic flow coefficient, and ∆Φ is the electrical potential difference between the two sides of the membrane. Before giving the explicit expressions for G and H for a finite-length cylindrical pore, we consider the extreme cases of a very thin membrane (L D) and a thick membrane (L D).
For L D, the membrane is regarded as a zero thickness sheet with a pore, where the entrance effect is dominant. The conductance is proportional to the diameter, namely G m = κD, where κ is the bulk electrical conductivity of the solution. If we restrict ourselves to the small Debye length regime (λ d D, as is true in the present MD simulations), the electro-osmotic flow coefficient is written as H m ∼ Dσ m λ d /µ, where µ is the solution viscosity and σ m is the charge density along the rim of the pore. For a long cylindrical pore (L D), the conductance is G c = πD 2 κ/(4L). For a small Debye length λ d D, the electro-osmotic flow coefficient is written as H c ∼ πD 2 σ c λ d /(4µL), where σ c is the surface charge density on the cylinder wall. Note that here and in what follows we use subscript m to represent quantities for the thin membrane with a pore (L D), and subscript c for the long cylindrical pore (L D).
We now give the expressions of G and H for a cylindrical pore of finite length L using the coefficients for the limiting cases mentioned above. The conductance G is simply obtained by combining the effects of the entrance and the cylinder in series: In the expression of the electro-osmotic flow coefficient H, we need to take into account the effect of internal pressure difference ∆ P, which is generated between the ends of the finite-length cylinder (see e.g. Ref. [89]). For convenience in the following discussion, we denote the electrical potential difference between the two ends of the cylinder by ∆φ. The internal values ∆ P and ∆φ are obtained using the two continuity equations for the current and flow rate. More specifically, the continuity of the current through the membrane is written as where the left-hand side of the second equality is the current outside the membrane and the right-hand side is that inside the membrane. The continuity of the flow rate Q z is written in terms of the flow rates generated by both the potential difference and the pressure difference [89]: where L m = −D 3 /(24µ) is Sampson's formula [90] for the permeance of a pore on an infinitely thin sheet and L c = −πD 4 /(128Lµ) is Poiseuille's law for the permeance of a cylindrical pipe. The left-hand side of the second equality is the sum of the flow rates outside the membrane induced by the electrical potential and pressure differences, and the right-hand side is that inside the cylinder. Now the explicit expressions of ∆φ and ∆ P are given by solving Eqs. (5) and (6), and the flow rate is then calculated using the first equality of Eq. (6) as follows: Effective pore diameter In our model of the electro-osmotic diode, the change in the current and flow is captured by introducing the effective pore diameter D eff , which replaces D in Eqs. (4) and (8) to capture the effect of the presence of a colloidal particle near the pore entrance. We propose the following form for D eff : where α > 0 is a model parameter (constant), P z (q, E z ) is a quantity proportional to the probability that the colloidal particle is close to membrane M 1 , and P xy (q) is related to the probability that the colloidal particle is around the pore; P 0 z and P 0 xy are the values at q = 0. Because P xy represents the colloidal particle motion in the x-y plane, it is assumed to be independent of E z .
In deriving this simple model expression, we assume the small variation in D eff is related to the probability p * (which represents the variables in the exponential in Eq. (9)) such that dp * ∼ −(D eff dD eff )/(πD 2 eff ), i.e. p * is in proportion to the ratio of the variation in pore area. Equation (9) is then obtained by integrating it under the requirement of D eff at p * → 0. Accordingly the essential feature of the diode is captured by this equation. If the colloidal particle is sufficiently far from M 1 and thus P z ∼ 0, then the effective diameter is equal to the pore diameter (D eff ∼ D). In contrast, if the colloidal particle is near M 1 (P z 0) and around the pore (P xy 0), then the pore is completely clogged (D eff ∼ 0).
In the following subsection, we discuss the analytical expressions for P z and P xy using the solutions of the Fokker-Planck equation.

Fokker-Planck equation for colloidal particle motion
Here, we discuss the motion of the colloidal particle employing the Fokker-Planck equations to obtain the expressions for the parameters P z and P xy , which are included in the effective diameter. The colloidal particle in the solution under an electrical field is subjected to the electro-phoretic force and random forces from the solvent particles, which leads to advection and diffusion.
In the following, the Fokker-Planck equations governing the probability density under this situation are separately formulated for motion along the z axis and in the x-y plane.
The Fokker-Planck equation for the probability density function p z is where the coefficient A(q, E z ) corresponds to advection and D is the diffusion coefficient for the particle, which we assume to be independent of q and E z . In the steadystate, the analytical solution is written as where z 1 and z 2 are the positions of boundaries at the surface of M 1 and M 2 , respectively (see dashed lines in Fig. 3).
Next, we focus on the motion of the colloidal particle in the x-y direction. We consider the existence probability distribution p r in the cylindrical coordinate system. The Fokker-Planck equation is written as follows: Here, we assume that coefficient A r is inversely proportional to r 2 , i.e., A r (q, r) = C A (q)/r 2 , to take into account radial advection due to electrical interaction between the charged particle and the pore. Then, the following analytical solution is obtained: The parameters in the effective diameter in Eq. (9) are defined as P z = p z (z 1 ) and P xy = δr 0 r p r (r )dr where δr is the radius of the vicinity region around the pore. We discuss the actual values of these parameters corresponding to the MD setup in Sec. with the aid of independent MD simulations conducted to estimate parameters.

Model parameters
Here, we determine the actual values of the model parameters for the system used for our MD simulations. We first estimate the diffusion coefficient D for the colloidal particle in the solution. To this end, we perform independent MD simulations for Brownian motion of the colloidal particle in the bulk electrolyte solution (see Supplemental Information S1 for details). From the mean square displacement of the particle, the diffusion coefficient is obtained as D = 1.7 × 10 −2 σ 2 0 /τ 0 . Because the colloidal particle charge has little effect on diffusive motion, we use this value throughout the following discussion.
We next consider the effect of advection in the z direction, which is incorporated in the model via the parameter A(q, E z ) in Eqs. (11) and (12). The motion of the colloidal particle with various values of charge q is tracked using an MD simulation setup where the colloidal particle is placed between two virtual membranes (same as M 2 used in Sec. ) with an electric field E z applied. (see Supplemental Information S2 for details). We obtain the probability density function for the colloidal particle along the z-axis, and determine the values of the coefficient A at each (q, E z ) by fitting the analytical solution given by Eqs. (11) and (12) using the value of D obtained above. We found the form A = βqE z with a constant β = −1.3 × 10 −2 σ 2 0 / 0 τ 0 . Finally, we examine the effect of advection in the x-y direction to determine the value of C A in Eqs. (14) and (15). Here, we consider the situation where the motion of the colloidal particle is constrained near membrane M 1 (see Supplemental Information S3 for details). This constraint is realized by placing the virtual membrane (same as M 2 ) at a distance of one diameter of the colloidal particle away from M 1 . To focus on the motion in the x-y plane, the pore charges are distributed as in Sec. , but the physical pore is not taken into account. Using the same value for D, we determine the value of C A by fitting the analytical solution given by Eqs. (14) and (15) to the existence probability density functions obtained with these constrained MD simulations. Accordingly, the value of C A was obtained in the form C A = γq, with γ = 6.2 × 10 −4 σ 3 0 /τ 0 q 0 .

Comparison with full MD simulation
With the obtained values of the model parameters D, A, and C A from independent MD simulations, we are able to estimate the effective diameter of the pore in the presence of the colloidal particle. The current and electroosmotic flow estimated by the model given in Eqs.  The other physical parameters in the model are set as follows. Parameters such as the surface charge densities σ c and σ m were determined from the geometrical setup. The electrical conductivity κ and viscosity µ, which are related to the properties of the electrolyte solution, could also be measured independently. In the comparison above, however, to bypass the independent measurements of these parameter values, we obtained them from a linear approximation of the current and flow rate in the region E z ≥ 0, where the pure current and pure electro-osmotic flow are approximately obtained because the colloidal particle is absent near the pore.

CONCLUDING REMARKS
In this study, we proposed an electro-osmotic diode that consists of colloidal particles between two membranes. The difference in pore density between two mem-branes makes the colloidal particles serve as nano-valves.
The rectification of the current and electro-osmotic flow is observed in MD simulations of the entire system. We established an analytical model that captures the key feature of the electro-osmotic diode, i.e., the colloidal nanovalve, where a colloid particle clogs the pore to prevent the backward current and flow.
The behavior of the colloidal nano-valve is modeled employing the introduced effective pore diameter, which depends on the existence probability density functions around the pore, as defined in Eq. (9). With the estimated model parameters, we obtained the forms of the probability density functions using the analytical expressions for the Fokker-Planck equations for the advection and diffusion of a particle. In our comparison with the full MD simulations, we estimated the parameter values using several independent MD simulations of colloidal motion, namely a bulk simulation to determine the diffusion coefficient, and two simulations for a colloidal particle in closed domains to determine the advection parameters in the z and x-y directions. Using the obtained parameter values, the full MD simulation results are successfully reproduced by the model equation, as shown in Figs. 2 and 5. This confirms that the proposed model reproduces the essential mechanism of the colloidal nano-valve. The quantification of the free parameter α is obtained via a top-down estimation in the comparison. Though one-time calibration of α is required for direct comparison, investigation of the diode performance for various situations becomes possible once the value of α is identified. Further investigation into the microscopic colloidal particle motion near the pore would allow us to determine the value of α using a bottom-up approach. This topic will be considered in future work. We also note here that this idea of effective pore diameter has other potential applications, such as molecular sieving, nano-filtration, and transports in porous media as in liquid electrolyte secondary batteries.
Throughout the paper, the membrane with a high number density of pores (membrane M 2 ) is regarded as a complete semi-permeable membrane, in which only the colloidal particle is affected by the potential barrier. This assumption was made to focus on the geometrical constraint of the colloidal particle. In a real membrane, however, the electro-osmotic flow is also driven in the pores in M 2 . This effect might not be negligible, and it would greatly enhance the magnitude of the current and the electro-osmotic flow while maintaining the basic mechanism of the nano-valves. Studies with the explicit configuration of M 2 would thus reveal this flow enhancement. Furthermore, the assumption of abstract materials for the colloidal particles, membranes, and electrolyte solution, can be replaced by the adoption of practical materials. Taking into account real materials, along with the design of appropriate experimental setups, is also a future research topic.
The proposed electro-osmotic diode is constructed with combining different physics, namely the electro-osmotic flow occurring in nano/micro pores, electrophoretic motion of colloidal particles, and the mechanical valve effect. We hope this crosscutting idea, including the mathematical modelling, will promote interactions among scientists in different fields. On the practical side, the present diode has potential application in rectifying flows under an AC electric field if the period is sufficiently longer than the characteristic relaxation time of the rectification. The nano-valves allow asymmetric flow and block backward flow almost completely if the physical parameters are appropriately chosen. Therefore, more efficient rectification compared to that achieved by existing AC electroosmotic systems is expected. We hope that the present simple configuration of the diode system will accelerate the development of practical electro-osmotic pumps in micro-and nano-fluidic devices. The bulk MD simulations are performed to investigate the diffusive motion of a colloidal particle. The colloidal particle is put in a cube simulation box as shown in Fig. 6. The system contains 6290 solution particles including 50 − q monovalent cations and 50 monovalent anions. The initial configuration is equilibrated such that the pressures for all directions are at 1 ± 0.1 0 /σ 3 0 . The colloidal particle with charge q has the same structure as that of Fig. 1 in the main text. No electric field is applied. The colloidal particle charge is set at q (= 6, 3, 2, 0 q 0 ), which are used in the full MD simulations in Fig. 2.
After long time has passed, the mean square displacement (MSD) r 2 = (1/3) x 2 becomes proportional to time t, which is a typical feature of diffusive motion. The MSD is then written in the form r 2 = 2Dt with D being the diffusion coefficient. The plot of the MSD for the cases of q = 6 q 0 is shown in Fig. 7. The diffusion coefficient D is obtained from linear fitting in the long-time regime. Practically, we used the data for t 20τ 0 , which is sufficiently larger than the mean free time. The colloid charge had little impact, and the mean value of the diffusion coefficient was obtained as D = 1.72 × 10 −2 ± 0.09 σ 2 0 /τ 0 . The simulations were performed for 10 7 steps with the time step dt = 0.002τ 0 . Three runs with different initial configurations were used for each q. In obtaining the MSD for each run, samples of for the time window of 100τ 0 (5 × 10 4 steps) were averaged with shifting the window by 1000 steps.

Advection of a colloidal particle along the electric field
The coefficient for the advection in the Fokker-Planck equation (A in the main text) is evaluated using the system in which both membranes M 1 and M 2 are the virtual membrane of potential barriers as shown in Fig. 8. The colloidal particle with the charge q has the same structure as in Fig. 1. The system contains 2410 solution particles including 20 − q monovalent cations and 20 monovalent anions, and the system size is L x = L y = 10 σ 0 , and L z ∼ 30 σ 0 . Here, as in the full MD simulation ions colloid electric field FIG. 8. The system of MD simulations. A colloidal particle is put between two potential barriers M1 and M2.
( Fig. 1), the system is equilibrated such that the pressure in the z direction P z is 1.0 ± 0.1 0 /σ 3 0 . The periodic boundary condition is assumed in all directions. The simulation system is illustrated in Fig. 8. Three production runs with different initial configurations are performed for q = 6, 3, 2 q 0 . At each production run, 10 7 time-step simulation is carried out with the time step dt = 0.002τ 0 .
From the trajectories of the MD simulations, the existence probability p z (z) of the colloidal particles are computed for the steady state. The applied electric field is varied in −0.7 ≤ E z ≤ 0.7 E 0 as in the main text. The distribution functions p z (z) are normalized such that z2 z1 p z (z )dz = 1, where z 1 and z 2 are the effective membrane surface, such that the positions of peaks of p z (z) coincide with the effective membrane surfaces.
The advection coefficient A is obtained using the numerically measured p z (z), by fitting the analytical solution of the Fokker-Planck equations Eqs. (11) and (12) with the least-squares method. The results are shown in Fig. 9. Here, the diffusion coefficient D is fixed at 1.72 × 10 −2 σ 2 0 /τ 0 , which was obtained in the Sec. . The obtained values of the advection coefficient A are plotted as a function of the electric field E z in Fig. 10(a). As shown by the fitted linear lines, the dependence on E z is linear at each value of q. Next, the slopes of the fitted linear lines are plotted as a function of q in Fig. 10(b), which is also linear. Altogether, we infer the form of function A(q, E z ) = βqE z , The coefficient is evaluated as β = −1.3 × 10 −2 σ 2 0 / 0 τ 0 from the slope of (b). Using this analytical form of A(q, E z ) together with the value of D, we can now obtain the analytical value of the existence probability at the surface of M 1 , p z (z 1 ), which is used in the definition of the effective pore diameter. In Fig. 11, we compare the analytically obtained p z (z 1 ) and the simulation results (symbols), showing a good agreement. 10  FIG. 9. The existence probability density function of the colloidal particle along the z axis for various values of the colloidal particle charge (a) q = 6 q0, (b) q = 3 q0 and (c) q = 2 q0. The symbols and lines represent MD simulation data and fitted curve, respectively. The electric field is varied from −0.7 E0 to 0.7 E0.

Motion of a colloidal particle in parallel to the membrane
Here, we discuss the existence probability of a colloidal particle close to the membrane. To this end, we consider a system with plane membrane as M 1 , and a virtual energy barrier M 2 that is placed near M 1 to restrict the motion of colloidal particles to the region close to System consists of a colloidal particle between physical membrane M1 and potential barrier M2 and solution. The interval between M1 and M2 is set to the extent of the colloidal particle diameter, such that particle motion is restrained in x − y plane.
the membrane. Specificaly, the position of M 2 is at z = L z /4 + e M1 /2 + 4σ 0 , where e M1 is the width of M 1 . In order to prevent the colloidal particle from being trapped in the pore, M 1 has no pore, while the pore charges are distributed in the same manner as in the main text. The numbers of solvent and ion particles are the same as in the simulation of Fig. 1. The system used here is shown in Fig. 12.
From the MD simulations we compute the distribution of colloidal particles in the r direction p r (r), where r is the radial coordinate in the cylindrical coordinate system in the plane parallel to the membrane. The tail of the probability distribution is appropriately corrected using the nature of the periodic boundary conditions. We plot the measured p r (r) in Fig. 13. As in the previous section, three production runs with different initial configurations are performed for q = 6, 3, 2 q 0 . At each production run, 10 7 time-step simulation is carried out with the time step dt = 0.002τ 0 .
We compare the results of MD simulation data with the analytical solution in Fig. 13. Note that the analytical solution is based on the negative point charge at the center of the membrane, while in the MD simulation the negative charges on the membrane are distributed in a finite region. Therefore, to pot weight on the tail (large r) in fitting the analytical functions, we consider the following objective function f : where h(r) is the values of the histogram obtained from the simulations, while p r is the analytical solution; δr ∼ 4σ 0 is a threshold, which is about the radius of the region where the direct effects of actual negative charges are significant in the simulation. The first term on the righthand side is the difference in integrals of the distributions for r < δr, and the second term is the difference at each point of r ≥ δr. Putting emphasis on the tail (r ≥ δr) of the distribution, we set a weight of w = 10. We used the downhill simplex method for the optimization. From the result of fitting, we found the dependence of C A on q to be linear, i.e., C A = γq. The value of the parameter γ was found to be γ = 6.2 × 10 −4 σ 3 0 /τ 0 q 0 . * h-yoshida@mosk.tytlabs.co.jp