Single-file transport of binary hard-sphere mixtures through periodic potentials

,


I. INTRODUCTION
Single-file dynamics refers to the collective motion in a many-particle system, where the particles cannot overtake each other because of their interactions and spatial confinements.Examples of such process are diffusion in pores of zeolites [1][2][3], colloidal particles confined to move to move in one-dimensional circular channels [4], particle motion through carbon nanotubes [5], in nanofluidic devices [6] and in mesoporous materials [7,8].Single-file transport often occurs in cell biology, e.g. in the movement of motor proteins [9,10], in protein synthesis by ribosomes [11,12], and in ion motion through narrow channels [13].
Recently, the Brownian asymmetric simple exclusion process (BASEP) was introduced as a paradigmatic model of single-file diffusion in a periodic potential [14][15][16].A striking feature of the BASEP is the sensitive dependence of its physics on the particle diameter σ.This is present even when ignoring a dependence of the particle mobility µ on σ.Such a dependence naturally occurs in fluid environments, where it can be accounted for by Stokes' friction law, µ ∝ 1/σ, for a single particle.More complicated forms may be relevant in the single-file dynamics of interacting particles due to hydrodynamic interactions.
Here we study the impact of polydispersity on the transport behavior in the BASEP.Polydispersity can hardly be fully avoided in experiments and it can be highly relevant in single-file transport, as, for example, in biological membranes and channels, in which transport of various types of proteins, ions, or other biomolecules takes place.Such membranes and channels are highly selective, with only a few types of particles being able to pass through.Underlying selection mechanisms often involve particle size.
We focus on a binary mixture of particles with hardcore interactions, i.e. two types of hard spheres with diameters σ A and σ B > σ A .Our primary interest is to understand how particle currents vary with the mixing ratio of the two types of particles and the hard-sphere diameters.Under single-file conditions, the particle currents also depend on the specific ordering of particles.We will show, however, that they exhibit a typical behavior for a given mixing ratio, and will therefore consider currents averaged over orderings.
We furthermore tackle commensurability effects, where the hard-sphere diameters are integer multiples of the period length λ of the external potential.Due to these effects, a basic unit cell can be defined in the σ A σ B -plane encompassing all hard-sphere diameters.By studying observables within this unit cell, one can predict their behavior anywhere else in the plane using transformation laws.
The transformation laws imply that it is interesting to study single-file transport of hard-sphere diameters under single-file conditions, where the diameters would actually allow the particle to overtake each other, i.e.where sums of two hard-sphere diameters are smaller than the pore diameter.By transformation to larger hard-sphere diameters, single file constraints become fulfilled.
The paper is organized as follows.After presenting the BASEP mixture model in Sec.II, we derive in Sec.III a mapping between probabilities of Brownian paths of many-body systems with particles of different hard-sphere diameters and discuss its implication on single-file constraints.The transformation between currents resulting from this mapping is provided in Sec.IV, which serves as a useful basis for explaining repeating features in the dependence of particle currents on particle size.Results from extensive Brownian dynamics simulations for the particle currents are presented in Sec.V. We interpret their dependence on the mixing ratio and particle diameters by considering the limit of a full coverage, where the currents can be calculated based on the Brownian motion of the center of mass in an effective potential.

II. BASEP MIXTURE MODEL
We consider a binary mixture of hard spheres A and B with particle diameters σ A and σ B ≥ σ A , which are driven through a narrow pore with a periodic structure in an aqueous environment as illustrated in Fig. 1.The overdamped Brownian motion of a sphere i of type α, α = A or B, is described by the Langevin equation where z i (t) is the position of the sphere, f α > 0 is a constant drag force, and U α (z) a periodic external potential reflecting the interactions with the pore walls; µ α and D α are the bare mobility and diffusion coefficient of particles of type α, where D α = k B T µ α .In accord with the Stokes law for the drag force, we set for the diffusion coefficients of the two particle types, where D is a bare diffusion coefficient of a reference sphere with diameter λ.The total number of particles is N , and the ξ i (t), i = 1, . . ., N , are Gaussian white noise processes satisfying ⟨ξ i (t For simplicity, we assume in this study the drag forces f α and the potentials U α (z) to be equal for both types of particles, Here, λ is the wavelength of the potential and U 0 the barrier between potential wells.The hard-sphere interaction implies the constraints For analyzing fundamental diagrams of particle currents as a function of particle densities, we investigate the particle motion in a closed system with M wells under periodic boundary conditions.Accordingly, the system length is FIG. 1. Schematic of a corrugated channel with a mixture of particles of two sizes (top).In the model, hard spheres of two diameters σA and σB are driven along a line by a constant drag force f through a periodic potential U (z) with wavelength λ (bottom).
Note that if the particle labels are assigned such as z 1 < . . .< z N at some time instant, the distance between the last and first particle must be smaller than The number density of particles is where N is the total number of particles, and N α and ρ α = N α /L are the numbers and number densities of particles of type α.The mixing ratio of the two particle types is Accordingly, ρ A = (1 − x)ρ and ρ B = xρ.Since the coverage N A σ A + N B σ B by the particles must be smaller than the system length L, we must have In the following, we adopt k B T , λ, and λ 2 /D as units of energy, length, and time, respectively.The potential barrier is set to U 0 = 6, i.e. it is six times larger than the thermal energy.
As the drag force f increases, the potential barriers gradually diminish, eventually vanishing at the critical force f c = πU 0 /λ, where the local maxima of the tilted potential [U (z) − f z] turn into saddle points.At f = f c , [U (z) − f z] has vanishing first and second order derivatives, which yields f c = πU 0 /λ.This critical force marks the transition from overcritical to undercritical tilting of the periodic potential.We are interested here in the undercritical regime f < f c and have set We simulated the time evolution according to Eq. ( 1) and subject to the constraints (4) by using Scala's eventdriven algorithm [17][18][19] with a time-step ∆t = 0.4 × 10 −3 .Particle currents J α were sampled after the system has evolved into a nonequilibrium steady state.All simulations were carried out for system sizes corresponding to M = 50 potential wells.

III. SINGLE-FILE CONSTRAINTS AND MAPPING BETWEEN PATH PROBABILITIES
The restriction of particle motion to one dimension in the model is motivated by considering the motion to take place in a narrow pore under single-file condition, where particles cannot overtake each other.If the pore diameter is d, single-file motion for both the A and B spheres would occur for 2σ A > d, as we defined σ B > σ A .Because σ B must be smaller than d, this implies the single-file condition Studying situations by the model, where this condition is violated, seems to be unphysical.However, as we will show now, it is still valuable to consider particle sizes that violate the single-file condition.This is because the current (and other quantities) in such systems can be mapped onto the current in a system with larger particle sizes that do satisfy the single-file conditions.Intuitively, such a mapping is imaginable due to the fact that a shifting of particle positions by a multiple integer of λ does not change the external forces in the periodic potential.Hence, one should obtain an equivalent system when increasing the distances between particles' centers of mass by integer multiples of λ.By increasing particle sizes, the interactions between particles can be kept the same, i.e. the new particle coordinates and sizes should obey the hard-sphere constraints (4).In addition, the change of positions and sizes should be accompanied by an appropriate change of the system length to accommodate all particles.
The mapping is strictly true when neglecting the difference between the diffusion coefficients of the two particle types as given by Eq. ( 2), and approximately when including this dependence, as we will see in Sec.V.
For D α = D, the probability P {z i (t)}; L, σ A , σ B of finding a path (trajectory) {z i (t)} in a system S with particle sizes σ α is equal to the probability ) This holds under the following conditions: (i) The path {z ′ i (t)} follows from the path {z i (t)} by a transformation specified below, see Eqs. ( 14) and ( 16).
(ii) The particle diameters differ by a multiple integer of the wavelength λ, where (iii) The particle numbers in both systems S and S ′ are the same, but the system lengths transform as This corresponds to a mapping between particle densities.
To show this, we proceed in two steps.First, we consider the case m A = m B = m.Second, we consider different m A and m B that are both even.By combining the results, we arrive at the conclusion in Eq. (10).

A. Case mA = mB
For m A = m B = m, we apply the coordinate transformation to the Langevin equations ( 1).The transformation leaves the external forces , and the hard-sphere constraints (4) become ( we obtain (z )/2, i.e. analogous constraints.
Moreover, the condition (z of the system length.Mathematically, the corresponding change of length ensures that the probability measure for the joint probability density of particle positions does not change. Thus, Eq. ( 10) is valid for m A = m B , if the particle diameters and system length are transformed according to Eqs. ( 11) and (12).

B. Case mA ̸ = mB
For m A ̸ = m B , we start by considering both m A and m B to be even.In that case, we apply the coordinate transformation r y 9 V k P x 9 y + P d S 7 z t H r V n u v t f t p p 9 n p V g N f w w u 8 x D a n u o 8 O P u C Q d U T 4 i k t c 4 Y e X e u f e h f f t L 9 W r V T n P c W N 5 3 / 8 A U H i d G A = = < / l a t e x i t > A < l a t e x i t s h a 1 _ b a s e 6 4 = " + y r P t m E q H 0 V w v x 9 o + F C p t o i i i 7 Corresponding regions in the σAσB-plane: The value of an observable at any point in the blue (orange) areas of the σAσB-plane can be inferred from the knowledge of the observable at a corresponding point in the blue (orange) area in the basic unit cell C bordered by the bold red line [see Eq. ( 17)].In the mapping between observables, the mixing ratio x is kept fixed, the density transforms according to Eq. ( 13), and the particle diameters according to Eqs. (11a), (11b) for σB > σA.Black circles illustrate a mapping for a point with σB > σA onto a point in C. A point with σB < σA is first reflected at the diagonal σB = σA, and is then mapped onto a point in C by applying Eqs.(11a), (11b).This is illustrated by the white circles, where the framed white circle indicates the mirror image. i.e.
Because all m i are even, this transformation again leaves the external forces f ext (z i ) = −dU (z i )/dz i unaltered.As for the hard-sphere constraints (4), we obtain ( λ of the system length, in agreement with Eq. ( 12).
Let us now consider a system with some diameters σ A < σ B .Applying the transformation for equal m α = m with m = int(σ A /λ) then yields an "equivalent system" ) is odd, a subsequent transformation with even m ′ A = 0 and even m ′ B +1 yields an unaltered σ ′ A ∈ [0, λ[ and a reduced σ ′′ B ∈ [λ, 2λ[.To conclude, the behavior of quantities in the σ A -σ B plane can be predicted from the knowledge in the basis unit cell More specifically, knowing the behavior of an observable in this basis unit cell of particle diameters for all mixing ratios x and particle densities ρ, one can infer the behavior of any observable at an arbitrary point in the σ A σ B -plane at some mixing ratio and density.Figure 2 illustrates the corresponding mapping.In the next section, we will demonstrate what this mapping implies for the particle currents, which constitute the primary focus in our study.
Turning back to the single file condition ( 9), particle diameters in the basis unit cell violating this condition are mapped to diameters in another unit cell in the σ A σ Bplane, where the conditions are satisfied.For example, by carrying out calculations for σ ′ A = 0.2λ and σ ′ B = 1.6λ, we can predict the behavior of quantities in a related system with particle diameters σ A = 2.2λ and σ B = 3.6λ satisfying σ B < 2σ A .

C. Mapping onto system of point particles
For particle diameters where σ A and σ B are both even or both odd multiples of λ, the mapping to the basic unit cell C yields σ ′ A = σ ′ B = 0, i.e. we obtain a system of point particles.
In a system of point particles, collective observables equal those in a system of independent (noninteracting) particles.Such observables are sums of equivalent terms over all particles, as it is the case for the total particle density and particle current.For example, for σ A = λ, σ B = 3λ, ρ = N/L, and any sequence with N α particles of type α, the current and the density profile is the same as in the system with N point particles (σ Hence, the quantities are also the same as in the system of noninteracting particles at density ρ ′ . Furthermore, the mapping to noninteracting particles allows us to gain a deeper understanding of the shape of the basic unit cell C in Eq. ( 17) (see also Fig. 2).In particular, one may wonder why σ B in the unit cell is not bounded by λ but by 2λ.
For hard spheres in one dimension, a mean interaction force between two particles in contact can be nonzero only if the mean external forces exerted on particle clusters left and right of the contact point are different, see [20] and the supplemental material of Ref. [16].If σ A and (σ A + σ B )/2 are multiple integers of λ, the external forces on all particles in contact are the same (for any position of the particles) and, accordingly, also the mean external forces exerted on clusters of neighboring particles in contact.This implies that the mean interaction forces vanish.The particles, however, keep their ordering.For obtaining the upper bound for σ B > σ A in the unit cell, one can determine the smallest σ B leading to zero mean interaction forces if σ A = 0.This is given by The argument can be generalized to three and more types of hard spheres with diameters σ 1 < σ 2 < . . .< The magnitudes of the currents and relative standard deviations are given by the colorbars.For completeness, we show results also for σB < σA (which are redundant in principle, as they mirror results for σA < σB).
The three σ A values are representative for small, intermediate and large diameters of the smaller A spheres in the basic unit cell.The high value of ρ = 0.8 is chosen to demonstrate the impact of the hard-sphere interactions on the behavior of particle currents upon mixing of the two particle types.We have also performed simulations for a small density ρ = 0.2 and found an analogous, but less pronounced variation of particle currents with x and the particle diameters.This finding let us conclude that the density is important for the magnitude of current variation with the mixing ratio but not for the overall functional dependence on x and the σ α .
Due to the single-file constraint, the currents are not the same for all particle orderings, and the question arises, how strongly currents fluctuate between orderings.To answer this question, we have calculated the relative standard deviation ⟨J 2 ⟩ − ⟨J⟩ 2 /⟨J⟩ (⟨J⟩ = J) as a function of x and σ B for parameters sets I and II, using 40 random sequences of N = 40 particles.The results in Figs.4(c) and (d) show that the relative standard deviation is not exceeding 10% for most values of x and σ B .Where it does, the current J = ⟨J⟩ is very low, as can be seen from Figs. 4(a) and (b), where we plotted J as a function of x and σ B for parameter sets I and II.If J is very low, even small changes of currents with the particle orderings can lead to large relative fluctuations.
To summarize, except for very low J, we can speak about a typical current for a given mixing ratio independent of the specific particle ordering.These typical currents, which we obtain by averaging over random particle orderings at fixed x will be discussed in more detail now.
Comparing the results for parameter sets I [D α = D, Fig. 4(a)] and II [D α = D/σ α , Fig. 4(b)], the location and overall pattern of regions with high and low currents in the xσ B -plane are similar.The main effect induced by the varying diffusion coefficients is to weaken differences in the current magnitudes.This demonstrates the minor impact of varying diffusion coefficients on current changes.
The mixing of A and B particles can lead to strong changes of the current by many orders of magnitude, note the logarithmic scale bars in Figs.4(a) and (b).Our identification of a basic unit cell C in Sec.IV is helpful to gain a first understanding of the patterning of regions with high and low currents.It implies a repetitive behav-ior of currents when σ B is changed by integer multiples of 2λ, and this is indeed reflected in Fig. 4. Namely, regions of maxima (dark orange) appear at σ B close to 2.4 and 4.4, and regions of minima of the current (green-blue) appear at σ B close to 1.7 and 3.7.Although there are differences in detail, the minima and maxima occur at similar positions for the three distinct σ A values in the basic unit cell.
The minima and maxima are most pronounced close to the full coverage, which is given by the maximal possible x = x max in the inequality (8), i.e. for σ B ≥ 1/ρ by where 1/ρ = 1.25 for ρ = 0.8 in Fig. 4. When decreasing x at fixed σ B in Fig. 4, the currents in general increase (decrease) if they have a minimum (maximum) close to x max .Hence, for understanding the overall variation of currents with x and σ B in the colormaps in Figs.4(a) and (b), we can focus on explaining the occurrence of the minima and maxima appearing near x max .
Let us first look at the local densities ρ α (z) = ⟨ i δ(z− z i (t))⟩/N α of particles of type α, where the sum runs over all particles of type α.We determined ρ α (z) from the simulations by calculating the mean number of α particles falling into a small bin around the position z.σ B ≃ 1.7 and 3.7?This can be done by considering that in the limit of full coverage, all particles in the system are in contact, i.e. they form a single particle cluster, and the current can be determined by the motion of this single cluster.For any given ordering of the particles, the equation of motion for the center of mass position z cm of the cluster is where and σ j,j+1 = (σ j + σ j+1 )/2.Accordingly, we can use z = z 1 for defining the cluster position.Its time evolution is given by the Langevin equation where μ = µ/N , D = D/N = k B T μ and ξ(t) is a Gaussian noise with ⟨ξ(t)⟩ = 0 and ⟨ξ(t)ξ(t ′ )⟩ = δ(t−t ′ ).For a particle ordering with diameters σ 1 , . . ., σ N the potential for the center of mass motion in Eq. ( 25) is where A = This potential is λ-periodic.The mean velocity v of the Brownian motion of a single particle in a periodic potential under a constant drag force is known exactly [21][22][23].
For the stationary current ρv we obtain (β = 1/k B T ) .
By averaging this current over particle orderings, we obtain J(x max , ρ, σ A , σ B ) = ⟨j(x max , ρ, σ 1 , . . ., σ N )⟩, ⟨. ..⟩ denoting the average over orderings.As an example, we show in Fig. 6 the so-calculated current at full coverage as a function of σ B for fixed ρ = 0.8 and σ A = 0.25 (solid line with crosses).In this numerical calculation, we have taken an average of 5 × 10 4 different particle orderings.The calculated values agree very well with simulation results that we obtained by averaging over 100 orderings.To check that the analytical results for full coverage are describing the limiting behavior of the FIG. 6. Mean particle current J(xmax, ρ, σA, σB) at full coverage as a function of σB for fixed ρ = 0.8 and σA = 0.25 (solid line with crosses) calculated from Eq. ( 27).An average over 5 × 10 4 different particle ordering has been carried out.The diamonds indicate simulations results for the same parameters ρ = 0.8 and σA = 0.25 at almost full coverage (see text).many-particle system, we performed these simulations by taking the parameters at full coverage and then adding one potential well to give the particles some free space to move.The corresponding simulation results indicated by the diamonds in Fig. 6 are well described by the analytical prediction for full coverage.Further insights into the behavior of J(x max , ρ, σ A , σ B ) can be gained by averaging the amplitude factor |A|, yielding the effective potential Note that we did not average u cm (z), as the different phases φ for the orderings would lead to interference effects.These phases, however, are irrelevant for the stationary currents.Figure 7(a) shows U eff (z) for the diameters σ B corresponding to minima and maxima of the current.For σ B = 1.7 and 3.7, the cluster has to surmount high potential barriers giving rise to low currents.By contrast, for σ B = 2.4 and 4.4 no potential barriers can be seen on the scale of the U eff axis in the figure, i.e. potential barriers are negligible compared to k B T , leading to high currents.
Closer inspection of U eff (z) allows us to identify the particle diameters leading to highest and lowest barrier heights at full coverage.For the density ρ = 0.8, these barrier heights

VI. SUMMARY AND CONCLUSION
We have studied Brownian transport of a binary mixture of hard-spheres A and B with two different diameters σ A and σ B ≥ σ A .The particles are driven through a periodic potential with wavelength λ under single-file conditions, where they cannot overtake each other.By extensive Brownian dynamics simulations we calculated particle currents for a wide range of mixing ratios x and particle diameters at several particle densities.A complex behavior in dependence of the mixing ratio and particle diameters was obtained, with currents varying by many orders of magnitude at high particle densities.Fluctuations of the currents between different particle reorderings turned out to be small compared to the mean values, except for the smallest particle currents.The strong variations of currents are only weakly affected when taking into account a scaling of mobilities with the particle diameter according to Stokes' friction law.
To understand the complex behavior, one can thus consider both particle types to have the same mobility.For that situation, we derived an exact mapping between probabilities of Brownian paths of different systems, where a system is related to another one with hardsphere diameters and system lengths differing by integer multiples of the potential wavelength.As a consequence, there exists a basic unit cell in the σ A σ B -plane, with σ A lying between zero and λ, and σ B between σ A and 2λ.Once the behavior of an observable is determined within the basic unit cell for all particles at a specific mixing ratio, it can be inferred for any combination of particle diameters and density at the same mixing ratio.
Application of the mapping to particle currents allowed us to understand that minima and maxima of them occur repeatedly in the σ A σ B -plane at σ B values separated by approximately two wavelengths.We further showed that the overall current variation with x and with the hardsphere diameters σ A and σ B can be inferred essentially from the current change with particle diameters in the limit of full coverage, when the sum of all particle diameters equals the system size.This change of current in the limit of full coverage can be explained from a calculation of an effective potential for particle transport.
The methods developed here for understanding Brownian dynamics of a binary hard-sphere mixture (mapping between Brownian paths, effective potential at full coverage) can be applied to arbitrary periodic potentials and extended to mixtures of more than two types of hard spheres.Following earlier studies for monodispersive systems [24], it will be interesting to explore how the collective transport properties are reflected in local transition dynamics of tagged particles.A further interesting extension is the investigation of systems, where the number of particles exceeds the number of potential wells.In such overcrowded systems, we expect that the particle transport is dominated by particle cluster waves, which at low temperatures appear as stable propagating solitons [16].
More elaborate analytical treatments may be possible in the future, as, for example, by approximate treatments of exact equations for the time evolution of particle densities and correlations derived from the many-body Smoluchowsky equation [14], by using dynamical density functional theory for hard-sphere mixtures [25,26], or power density functional theory [27,28].Based on earlier investigations [29,30], we believe that an understanding of hard-sphere mixture dynamics in periodic structures will provide a suitable basis for treating systems with other particle interactions.
Our theoretical treatment of driven Brownian motion of a binary mixture could lay a suitable basis also to understand anomalous uphill diffusion against a concentration gradient of one particle type, while the other type of particle behaves normally.This effect was experimentally observed in two-component diffusion through nanopores [31] and was related to cross-coupling terms (nondiagonal Onsager coefficients), i.e. that a concentration gradient in one component influences the current of the other component.We believe that nondiagonal Onsager coefficients are strongly affected by changes of the particle size.

FIG. 4 .
FIG. 4. Stationary currents [panels (a), (b)] and their relative standard deviation for different orderings of particles [panels (c), (d)] in dependence of the mixing ratio x and particle diameter σB for the parameters sets I [Dα = D, panels (a),(c)] and II [Dα = D/σα, panels (b), (d)] (ρ = 0.8, three representative values of σA).The magnitudes of the currents and relative standard deviations are given by the colorbars.For completeness, we show results also for σB < σA (which are redundant in principle, as they mirror results for σA < σB).

Figure 5
FIG. 5. Averaged density profiles for tight particle packings.In (a) particles displace each other greatly, while in (b), the particles are confined in to potential wells blocking each other from motion.Gray dashed line indicates potential minimum.
β(u cm (y)−N f y−u cm (z)+N f z)] are shown in Fig. 7(b) as a function of σ B for various σ A .Low and high barriers occur at similar σ B for different σ A , which explains the aforementioned insensitivity of the locations of regions with low and high current in the σ A σ B -plane [see the discussion of Figs.4(a) and (b)].

FIG. 7 .
FIG. 7. (a) Effective potential U eff (z) for the single cluster formed by all N particles in the limit of full coverage at values of σB corresponding to minima and maxima of the current in Figs.4(a) and (b).(b) Heights of the potential barrier to be surmounted by the cluster as a function of the particle diameter σB for various σA.