{'fraction_non_alphanumeric': 0.029736286094373583, 'fraction_numerical': 0.008607872290476564, 'mean_word_length': 4.701204550524203, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of Fluid Mixtures in Nanospaces

A multicomponent extension of our recent theory of simple fluids [ U.M.B. Marconi and S. Melchionna, Journal of Chemical Physics, 131, 014105 (2009) ] is proposed to describe miscible and immiscible liquid mixtures under inhomogeneous, non steady conditions typical of confined fluid flows. We first derive from a microscopic level the evolution equations of the phase space distribution function of each component in terms of a set of self consistent fields, representing both body forces and viscous forces (forces dependent on the density distributions in the fluid and on the velocity distributions). Secondly, we solve numerically the resulting governing equations by means of the Lattice Boltzmann method whose implementation contains novel features with respect to existing approaches. Our model incorporates hydrodynamic flow, diffusion, surface tension, and the possibility for global and local viscosity variations. We validate our model by studying the bulk viscosity dependence of the mixture on concentration, packing fraction and size ratio. Finally we consider inhomogeneous systems and study the dynamics of mixtures in slits of molecular thickness and relate structural and flow properties.


I. INTRODUCTION
In recent years, the experimental study of the flow of liquids in nanometric structures has become an important and rapidly evolving discipline driven by technological applications, such as lab-on-a-chip systems, electrophoresis and electro-osmotic pumping. The mixing and separation of different molecules according to their size or chemical and transport properties is a major issue in chromatography, nanofluidic logic gates and electrokinetics. Many applications are foreseen also in biotechnologies, because very small quantities of liquid components are sufficient for analysis and synthesis, including the possibility of isolating and analyzing or manipulating a single molecule in a single nanochannel or surface with nanoscale features [1,2].
These technological advances require a better theoretical understanding of the fluid properties in confined geometries and under non-equilibrium conditions. At the nanoscale, new areas of physics, chemistry and materials science come in. The systems become more surface-like and molecules never explore a bulk-like environment and because inhomogeneities have a great impact some interesting phenomena arise. Whereas a large body of information has been accumulated in the last thirty years concerning the physics of inhomogeneous fluids under thermodynamic equilibrium conditions [3], the state of the art of flowing fluids is not so advanced [4][5][6][7].
It is crucial to understand how the interplay between the various forces occurs and how it changes when going from macrosystems to nanosystems. It is clear from dimensional considerations that the relative importance of forces changes with the typical dimensions of the systems and we expect that the smaller the sample under scrutiny the more important the role of surface forces is. Not only nanoscale flows are dominated by viscous effects if turbulence is absent, but wherever frequent molecule-channel wall collisions occur besides molecule-molecule collisions, we expect a modified frictional resistance [8].
The generalization of fluid dynamics from pure to multicomponent fluid mixtures composed of different components or species requires the introduction of some new concepts. While traditional computational fluid methods involve some ad hoc extrapolation of Navier-Stokes equation to confining geometries, our approach is based on a microscopic kinetic equation incorporating the effects of the inhomogeneities in its structure.
The major difference with respect to one-component fluids is the phenomenon of concentration diffusion. We shall study a multicomponent kinetic equation [9][10][11][12][13] that has been considered in the past by several authors with the specific goal of computing the transport coefficients from a microscopic approach. Many of these studies, which were based on the Boltzmann equation or on its extension to the dense case, the Boltzmann-Enskog equation, require some analytical and numerical effort, so that a simpler treatment of the interactions has been proposed. A very popular approximation is represented by the phenomenological Bhatnagar-Gross-Krook (BGK) equation whose merit is to reduce the complexity of the original Boltzmann collision kernel Ω αβ between species α and β by introducing a simple relaxation-time ansatz [14][15][16][17][18]. The BGK, in spite of being even less accurate than the original Boltzmann equation, has enjoyed a great popularity especially in conjunction with the Lattice Boltzmann method (LBM) thanks to the simplicity of the collision kernel leading to a considerable speed-up of the numerics [19] .
However, in the multicomponent case the choice of the form of the BGK relaxation term is not unique and hard to infer from the original Ω αβ , and the literature abounds of several proposals [20][21][22][23][24].
In order to describe with sufficient accuracy the fluid structure at length scales comparable with the size of the particles we shall resort to methods similar to those of density functional theory (DFT) employed in the study of equilibrium and non equilibrium properties [25][26][27][28]. In the case of hard-core fluids, DFT and its dynamical extension gives excellent results and can be extended to more realistic fluids by using the van der Waals picture of decomposing the total inter-particle potential into a short-range repulsive potential and a long-range attractive potential tail. The first is treated by means of a reference hard-sphere system whilst the second is considered within the random phase approximation (RPA). Such a decomposition can describe phase separating systems, wetting phenomena, two phase interfaces, but predicts transport coefficients of purely hard-core systems [29].
The hard-core part of the interaction is treated within the revised Enskog theory (RET) [30] that considers the non-local character of the momentum and energy transfer and takes into account the static spatial correlations among particles. This level of theory does not incorporate the time correlations which are responsible for memory effects and hydrodynamic contributions to the transport coefficients [31].
With respect to the dynamical DFT [32,33], the present theory preserves Galilei invariance and therefore displays full hydrodynamic behavior in the limit of slowly varying fluctuations about the equilibrium reference state.
With the aim of deriving a practical numerical approximation to study mixtures in inhomogeneous situations, we extended the method of Dufty and coworkers [34] to the multicomponent case. The method is a compromise between the RET, of which it retains the accuracy as far as the momentum and energy transfer are involved, and the much simpler BGK, which we employ to evolve the non-hydrodynamic moments of the distribution functions [35][36][37].
The final product of our theory is a coupled system of simplified equations for the density distributions of individual species describing both the streaming and the collisional stages. These equations are solved by an appropriate extension of the LBM algorithm to include both particle-particle interactions and particle-wall interactions. A simple analysis of the equations is used to derive explicit expressions both for equilibrium thermodynamic quantities, such as pressure, compressibility, etc., and for non equilibrium transport coefficients. It is important to stress the difference between the present work and other LBM based approaches for non-ideal fluids, where only attractive interactions are accounted for by the so-called Shan-Chen pseudo-potential, whose justification is purely mesoscopic, while transport coefficients enter as free parameters of the theory [38,39].
Interestingly, over two decades ago Davis and coworkers proposed on phenomenological grounds a local average density model (LADM) of viscosity and diffusivity, based on the principle of local averaging of the density over the volume of a hard sphere [40]. The present approach derives similar expressions for the transport coefficients, but differs from the LADM in the way it solves the resulting evolution equations. Pozhar gave a more rigorous foundation to the transport theory in non-uniform systems, but the resulting theory could not be developed into a computationally tractable method [41,42]. In the present work, we claim that our microscopically founded theory lends itself to relatively simple numerical solutions based on the use of the LBM even in the presence of complex geometries. The numerical application regards fluids confined to a spacing of a few monolayers between solid surfaces, where they may behave very differently from bulk conditions and exhibit unusual properties such as a dramatic increase in viscosity.
The paper is organized as follows: in Sec. II we introduce the model and the evolution equations for the individual distribution functions of the multicomponent system. and in Sec. III we discuss the approximations involved in the theory. In Sec. IV we test the predictions of the resulting evolution equations using a simple confining geometry. In Sec. V we validate the resulting equations numerically and simulate the flow of a mixture in narrow channels. Finally in Sec. VI we present our conclusions and perspectives. Two appendices are added to the paper in order to illustrate the relationship between the collision kernel and the non-ideal part of the chemical potential and discuss some technical details related to the discretization procedure.

II. KINETIC EQUATIONS FOR THE DISTRIBUTION FUNCTIONS
We consider a fluid mixture of N α particles of mass m α , with α = 1, M , where M represents the number of components, mutually interacting with the potentials U αβ (r, r ′ ). The probability density of finding a particle of type α at position r with velocity v at time t is proportional to the the single-particle distribution functions, f α (r, v, t), normalized as follows: (1) The evolution of f α is given, in principle, by the exact BBGKY hierarchy of dynamical equations for the reduced distribution functions [9], whose first member is the following kinetic equation for the α component: where F α (r) is an external velocity independent force field acting on component α and Ω αβ represents the effect of the interactions among the fluid particles of type α and β. In principle, the exact form of Ω αβ is known, but it involves the two-particle distribution functions, which in turn depend on the three particle distribution functions. In order to close this hierarchy different approximations have been proposed, the first historically being the Boltzmann equation expressing Ω αβ in terms of products of single-particle distributions. Other specific choices of Ω αβ in eq.(2) correspond to different kinds of approximations and lead to equations such as the BGK equation [14], the Vlasov or RPA equation [43] and the RET equation [30].
We seek a description in terms of a small number of macroscopic variables which are functions of the fixed position r and the time t. The procedure leading from eq.(2) to the evolution equations in terms of these variables, closely follows the analysis of the one-component case [36] and allows to derive a set of hydrodynamic equations, in the case of a M-component fluid mixture. In D spatial dimensions, there exist L = (D + M + 1) conservation laws, for the mass of each species, global momentum and energy. Thus, one can consider the evolution of L linear combinations of velocity moments of f α which are identified with the hydrodynamic fields of the mixture.
In the following, we shall adopt the Einstein summation convention for repeated Cartesian indexes, while sums over components are written explicitly. Let us first write the statistical expressions for the mass density of component α: where n α is the associated number density. The total average mass density is given by: while the total number density is n ≡ α n α . The first velocity moment of the phase-space distribution corresponds to the average local velocity, u α , of the α component We also consider the local barycentric velocity: and the so called dissipative diffusion current, measuring the drift component α with respect to the center of mass velocity: and having the property α J α (r, t) = 0.
B. Balance equations.
The macroscopic variables ρ α , u α , u obey a set of evolution equations derived from eq. (2) as follows. Upon multiplying by m α and integrating over the velocity degrees of freedom both sides of eq. (2) one obtains the continuity equation for the mass density of the species α: After summing over components, one obtains the global mass continuity equation: In establishing eq. (8) we have used the property that the interaction does not affect the partial number densities, i.e. dvΩ αβ = 0. The momenta of the single species are not locally conserved quantities, due to the interaction between different components. To derive the associated governing equations we first multiply eq. (2) by m α (v − u) and then integrate over the velocity with the result: In analogy with pure fluids, in eq. (10) we have introduced the kinetic contribution to the partial stress tensor: Notice that in order to ensure that π α ij is independent of a superimposed uniform translation of the entire flowing fluid its definition contains the difference (v α i − u i ). In eq. (10) we have also defined the projection of the collision term onto the velocity: representing the contribution to the stress tensor arising from interactions. This integral will be explicitly computed in Sec. III A for a specific model. Finally, the term ρ α w α i w α j in the l.h.s. of eq. (10), has nothing to do with viscous effects, but represents an additional stress of dynamical origin due to the different average velocities of the two components [44] and vanishes in the pure fluid. Summing over species eq. (10) we obtain the local conservation law for the total momentum, ρu, with the total kinetic component of the stress tensor given by: The set of eqs. (8) and (13) is not determined unless additional relations specifying the pressure tensor (14) in terms of the selected macroscopic fields are given. Hydrodynamics at the Navier-Stokes level closes the equations by assuming a simple linear relation between fluxes and derivatives of the fields via phenomenological transport coefficients (such as the viscosity and conductivity) whose values are not known a priori, but which must be taken from experiments. Kinetic theory, working at the level of the non-equilibrium distribution functions, avoids such assumptions.

III. SIMPLIFIED KINETIC MODEL
In a nutshell, our strategy consists in solving the kinetic equations rather than the balance equations (8) and (13). It involves concepts borrowed from the lattice Boltzmann method (LBM) and from density functional theory (DFT). The LBM can be regarded as an efficient approach to treat inhomogeneous systems and multiphase flows and provides an alternative to the solution of the Navier-Stokes equations. Conditions at solid boundaries are relatively easily implemented using the LBM, making it well suited for handling complex geometries. The LBM becomes particularly useful to treat non-local effects, such as wetting layers, interfaces, layering etc., if used in conjunction with a simplified kinetic approach suggested by Dufty and coworkers in the case of one-component fluids [34].
To introduce the simplified kinetic model for mixtures, we refer to ref. [36], where we discussed the simpler onecomponent case and reduced the full transport equation (2) to a simpler form. The extension to the multicomponent case requires some subtleties which are needed in order to satisfy a key requirement of this approximation: the modified kinetic equation and the original one must share the same equations for the partial densities and momenta, i.e. the distributions must have the same low order projections. Doing so, we are able to treat exactly the contribution of the collision operator to the hydrodynamic equations. As a result, we need to make approximations only on that part of the collision operator which acts on the higher order, non-conserved modes. This latter contribution can be handled in the spirit of some cruder approximation, such as the BGK-like approximation [18]. Since the velocity distributions in the final equilibrium state are known it is easy to construct the BGK collision operator in such a way that it verifies some basic properties of the original problem, namely particle and momentum conservation. Our two-component extension of the simplified kinetic equation reads: where ω is a collision frequency which depends on the two species and and The proposed form (15) differs from the existing literature about mixtures in a few aspects. In the first place, the mixtures were treated either by including the full collision operator or by assuming the BGK to be a valid replacement. The simplified kinetic method, to the best of our knowledge, was never applied to mixtures. Our choice for the BGK term is based on a single relaxation time for the two components and represents perhaps the minimal model satisfying the conservation laws compatible with the Dufty and coworkers strategy.
A few comments are in order. First of all, eqs. (8), (10) and (13) are immediately recovered after multiplying eq. (15) by m α , m α (v−u) and integrating over the velocity, v. Second, the coupling between the evolution equations of the two distributions f α occurs through the nonlinear collision kernel and this fact is also reflected in the BGK representation of that kernel. Equation (2) contains two terms Ω αβ , one for each type of interaction, but in the modeling of the BGK relaxation time approximation it is convenient to choose a single collision term. This assumption is valid not only because the coupling between the two distributions is achieved via the barycentric velocity u contained in ψ α , but also through the term C αβ . Third, the classical BGK approximation when applied to a one-component fluid would involve the difference between the distribution f α and the local Maxwellian ψ α . In the present case a non exponential prefactor in eq. (17) has been introduced and serves as to "orthogonalize" the term ω[f α −ψ α ⊥ ] to the "true" collisional terms, C αβ , so that the BGK does not lead to modifications of the hydrodynamic equations. Also notice that the first velocity moment of ψ α is n α u, whereas the first moment corresponding to ψ α ⊥ is n α u α . The present ansatz reduces to our one-component equation [36] in the limit ψ α ⊥ → ψ α . The above choice satisfies the indifferentiability principle which dictates that when all physical properties of the species are identical, the total distribution f = f A + f B obeys the single species transport equation.

A. Interactions
We consider the interaction term featuring in eq. (2) which according to the BBGKY hierarchy can be written as where f αβ 2 (r, r ′ , v, v ′ , t) represents the two-particle distribution function and the U αβ ((r − r ′ ) are the pair potentials between particles of species α and β.
In general, the intermolecular potential, U αβ , contains both a repulsive and an attractive contribution. The latter plays a major role in determining the thermodynamic properties of the system. In the following we shall separate the repulsive and the attractive contributions by representing the short range repulsion between the particles through repulsive hard sphere potentials, characterized by the three diameters σ AA , σ BB and σ AB = (σ AA + σ BB )/2 plus an attractive contribution: The repulsive part is treated within the following RET approximation [45]: where The quantities g αβ (r, r ± σ αβ k) are the three hard sphere pair correlation functions evaluated when the particles of species α and β are at contact. The RET collisional term describes the large momentum binary exchanges due to close encounters between the cores, whereas the attractive term can be seen as a slowly varying average force, involving small momentum exchanges. In detail, the attractive tails of the potentials are treated in a mean field approximation, rewriting the second term of the r.h.s. of eq.(19) by factorizing the two-particle distribution functions that is, by neglecting velocity correlations completely: where g αβ (r, r ′ , t) is the inhomogeneous the pair correlation function of the reference hard-sphere system. Integrating w.r.t. the velocity, we obtain the RPA approximation for the attractive term: where we introduced the molecular fields In order to evaluate the interaction contribution to the momentum balance equation, we use the definition (12) together with eqs. (20) and (24). After some algebra, using explicitly the collision rules eqs. (21) we obtain the following expression: where we have introduced the reduced mass µ αβ = m α m β m α +m β . To proceed further, we invoke a local equilibrium approximation and replace the distribution functions within the integrals (26) by local Maxwellians and perform the integrals involved in eq. (20), by expanding u α (r, t) about u(r, t). We also expand the functions at point (r + σ αβ ) about their values at r up to first order in the differences. The procedure is completely analogous to the one followed in ref. [36] and leads to the final result: B. Decomposition of the interaction term into specific force contributions.
It is instructive to rewrite the term featuring in formula (28) as a sum of more specific forces: We identify the force, F α,mf , acting on the α particle at r due to the influence of all remaining particles, being the gradient of the so-called potential of mean force: a drag force with local character: and a viscous force of non local character: We have further defined a friction tensor: and a viscosity kernel: whose relation with the transport coefficients will become clear in the following. The term F α,mf is related to the gradient of the excess chemical potential of species α, µ α int (r), due to the interactions among the particles. In Appendix A we show that in the case of slowly varying density profiles one has the result The drag term is proportional to the velocity difference of the two species and is present even in the absence of velocity gradients. The viscous term, on the other hand, is proportional to velocity gradients. We can consider the case u A = const and u B = const and u A = u B and uniform densities so that the tensor γ αβ ij becomes isotropic By using the definition (33) we obtain It is interesting to show how the microscopic viscosity kernel (34) is related to the viscosity coefficient. Let us assume a simple case u A = u B = u with u(r, t) = (0, 0, u z (x, y)) (38) and uniform densities, n A and n B . In this case the only non vanishing component of the viscous force is directed along z. After some lengthy, but straightforward algebra (see ref [36]) we obtain: which allows to identify the collisional contribution to the shear viscosity coefficient by using the macroscopic relation from which we can write the explicit formula: In the case of pure fluids, this expression of the viscosity reduces to the one predicted by Longuet-Higgins [46], which is known to describe successfully fluids of hard molecules even when the radial distribution function is momentarily isotropic. In spite of the fact that the trial distributions correspond the the local equilibrium state, the singular nature of the interactions allows a finite flux of momentum. Notice that the RPA attractive term does not contribute to the transport coefficients, but only to the potential of mean force.
C. Fischer-Methfessel prescription for the pair correlation functions of a non uniform system.
We have deduced the formulas for the viscosity and drag coefficient in the case of systems of uniform densities. In order to apply our theory to non-uniform cases we need to specify how the pair correlation functions are modified with respect to the bulk. The pair distribution function g αβ appearing in the above formulas is constructed according to a generalization of the prescription put forward by Fischer and Methfessel (FM) [47]. One first defines the smeared densitiesn α (r) via a uniform averaging of n α (r) over spheres of volume V α = πσ 3 αα /6 centered at the point R αβ located between the centers of the two spheres at r α and r β , respectively: The smeared densities are [48]:n The contact values of the pair correlation functions are computed using the bulk expressions obtained from an extension [49,50] of the Carnahan-Starling equation to mixtures: where the functions ξ n are evaluated for the values of the smeared densities at point R αβ between the two particles at contact:ξ The pair correlations at contact are evaluated in the inhomogeneous system using the generalized FM ansatz: From the knowledge of g αβ , one can immediately obtain the equation of state for the mixture: After having made the connection between the collisional term C αβ (r, t) and the three types of forces acting in the fluid eq. (29), we can use the large body of knowledge accumulated over the last twenty years about entropic forces in multicomponent systems [51]. These forces, for instance, determine an effective attraction towards a flat wall of the particles with the larger radius induced by the presence of the small particles. For similar entropic reasons two large particles immersed in a sea of small particles experience an effective attraction. To the best of our knowledge, the interplay between entropic forces and viscous forces is not very well known and is worth to be explored. As an order of magnitude estimate entropic forced are typically k B T nσ 2 AB , whereas the ratio between viscous forces and entropic forces is Although recent versions of equilibrium DFT are more accurate than the FM method, it is more convenient to work within the FM approximation for a series of reasons: a) the FM approximation is very simple as compared to the Rosenfeld [52] prescription for pair correlations; b) our theory not only requires the effective potential ∇µ α int , but also the viscous and frictional forces which are not given by the DFT method; c) computationally the FM method is much faster, although it becomes numerically unstable at high packing fractions.

D. Numerical validation
We consider a binary mixture of HS with σ AA = 4, σ BB = 8 and composition X A = n A /(n A + n B ). The first validation of the numerical algorithm, succinctly described in Appendix B, is obtained by calculating the shear viscosity of the binary system and initially subjected to a sinusoidal shear wave of small amplitude in bulk conditions. The measured decay time is found to have the same temporal decay for both species, and the resulting viscosity curves are reported in Fig. 1. The data are compared with the analytical expression for the decay time obtained by solving the linearized hydrodynamic equations and using the theoretical viscosity of eq. (41). Notice that the viscosity increases monotonically with the packing fraction ξ 3 of the system and depends on the molar fraction X A . At fixed packing fraction, the larger the concentration of the big particles the larger is the viscosity.
The agreement between numerical and analytical data is remarkable. At large packing fraction, a degree of departure from the analytical data indicates the major role played by the collisional integral that apparently is not fully captured by the numerical quadratures. This is an effect of the discretization and of the linear interpolation employed in order to compute some off-lattice points used in the calculation of the surface integrals. From the viewpoint of numerical stability, the range of packing fraction handled in simulation shows that the method is stable up to packing fractions of ξ 3 < 0.4, rendering the approach usable in many experimental condensed matter contexts. In proximity of hard walls and under strong flow conditions, the numerical stability range reduces to packing fraction below ξ 3 ∼ 0.3.

IV. FLOW IN A SLIT-LIKE PORE. TOY MODEL.
We now consider a binary mixture of particles of equal masses and different sizes (σ AA < σ BB ) confined in a narrow slit-like pore represented by two parallel smooth plates having infinite area and separated by a distance H + σ AA (see fig. 2). The walls are parallel to the yz plane and a flow along the z direction is induced by the presence of a uniform external field, F, also directed along z. Since H is the typical size of the system, u the average flow velocity and η the viscosity, the Reynolds number is Re = ρuH/η. In the small Re regime, flows have negligible inertia forces. The viscous and pressure forces should be in approximate balance and a reduction, similar to the one occurring between the macroscopic Navier-Stokes and the Stokes equation, takes place also at the present microscopic level of description.
The first remark is that, besides being small in systems of microscopic size and in the limit of low flow velocities, the non linear term in the momentum equation vanishes. In fact, if one imposes the boundary conditions typical of a Poiseuillle flow, in a straight channel, it vanishes due to symmetry.
We rewrite the l.h.s. of eq. (10) as: The first parenthesis in the r.h.s. is zero due to the continuity equation of species α, the second parenthesis is also zero due to the geometry of the problem and the stationarity of the flux, and the last term vanishes because w is directed along z and ρ α varies only along x. The flow induced by the presence of uniform fields parallel to the walls is characterized by velocities along the z direction u α = (0, 0, u α z (x)) obeying the following set of stationary equations for each species: In order to gain insight into the problem, we first set up and solve a toy model, obtained from the full model under some crude approximation. We replace the true non-uniform density profiles by two slabs of constant densities, n A , n B , respectively: As shown in ref. [36] the kinetic contribution to the viscosity due to the individual species is η (K) AA = kB T ω n A , while the collisional part is given by formula (41) . Notice that the above prescription for the kinetic contribution to the viscosity is a consequence of our single time relaxation ansatz. Let us introduce the abbreviations: When the two species have different diameters, the smaller species can approach the wall up to a distance σ AA /2, whereas the larger B species can only reach the distance σ AA /2 + a, where a = (σ BB − σ AA )/2. In the following, we take H > 2a and measure distances, x, from the position of closest approach of particles A to the wall as illustrated in Fig. 2.
We now determine the velocity profiles consistent with our toy model. In the central region a < x < (H − a) we obtain the following coupled differential equations for the partial velocities: In the excluded layer for the big spheres, that is for x < a or x > (H − a), the equations read: The solution with boundary conditions of zero velocity u A z at x = 0 and x = H, together with zero velocity u B z at x = a and x = H − a, and continuous u A z at x = a, reads: while for a < x < (H − a) is given by: and with and The main result of this illustrative model is the presence of different velocity profiles for the two species. The velocity profile of the bigger particles is generally lower than the corresponding profile of the smaller particles, because B effectively sees a narrower slit, so that the confinement effect resulting in a Poiseuille parabolic-like profile is more relevant than for the other species. The velocity difference between species remains finite in the center of the pore. Fig. 3 displays the velocity profiles of the two components for two different values of the composition and at low packing fraction. We will see below that even such a simplified model is able to capture some important features which are better studied by means of our full numerical method.

V. LBM NUMERICAL RESULTS
In the present section we solve numerically the equations for the phase space distributions and study the associated steady state properties. The most interesting applications of our method are those which regard truly inhomogeneous situations where one observes the interplay between the microscopic structure and the flow properties. The small amount of material in the confined region makes experimental probing of confined fluid extremely difficult and thus numerical methods are welcome.
We reconsider the Poiseuille problem of the previous section using the full power of our theory to obtain the density and velocity profiles of the two components when a uniform field parallel to the plates is applied.
The available region for each hard sphere can be computed according to the density profiles that show where the centers of the spheres can range. Given that the sphere radii are multiples of the mesh spacing, the extrema of the available region can only fall on a mesh point. In this respect, the values of H are equal to 7, 13 and 29 mesh-points for the three channels, that is H/σ AA = 1.75, 3.25 and 7.25 respectively. On the other hand, for the velocity profiles, we have added a mesh spacing on the left and right extrema to account for the non-slip boundary condition as imposed via the bounce-back scheme in the numerical method. When looking at the intervals from the velocity point of view, they look wider. Fig. 4 displays results for the density and velocity profiles for each species in a channel of width H = 1.75σ AA for a mixture of particles with diameters σ AA = 4 and σ BB = 8 in two different cases: packing fraction ξ 3 = 0.084 and composition X A = 0.25 and packing fraction ξ 3 = 0.073 and composition x A = 0.75. In the upper left panel the mixture is almost a gas and the density profiles for both compositions considered do not display any significant structure, whereas in the lower left panel the velocity profiles have nearly parabolic shapes, each species having its own curvature. One may conclude that in such a regime the flow behavior of the confined fluid is Poiseuille-like. As we increase the packing fraction to ξ 3 = 0.42 (composition X A = 0.25) and packing fraction ξ 3 = 0.26 (composition X A = 0.75), we observe the behavior reported in the right panels of Fig. 4. The large species B is unable to display significant peaks, due to the small plate separation that inhibits the formation of two layers of B particles. The small species, instead, starts developing some density enhancement in the proximity of the walls (upper right panel). The velocity profiles instead still bear a parabolic behavior with the smaller species always being faster than the bigger species.
In Fig. 5 we consider the same state conditions but a wider slit, H = 3.25σ AA , such that more than one layer of particles B can fill the gap. In the low packing region, at ξ 3 = 0.073 and at ξ 3 = 0.084 (upper left panel), both density profiles are rather flat and the associated velocity profiles do not display significant structure, the smaller particles being faster (lower left panel). As we consider larger packing fractions, ξ 3 = 0.42 and composition X A = 0.25, we observe a stronger structure in the large particles than in the small particles. The other packing fraction, ξ 3 = 0.26 (composition X A = 0.75), also produces some significant peak structure at the walls for the large particles, but the small particles have a rather flat density profile. The small HS are somehow enslaved by the large component as far as the structure is concerned. In the velocity profiles (lower right panel) we observe that the B component for X A = 0.25 has a higher velocity than the small counterpart near the center-line of the slit. This effect is determined by packing: small particles have more room near the walls while large particles move faster near the center-line.
In Fig. 6 we consider a channel with H = 7.25σ AA . In the low packing region (left panels), there is not significant structure, the velocity profiles are Poiseuille-like and display the same trend as in the corresponding panel of Figs. 4 and 5. More interesting is the situation at higher packing, shown in the upper right panel. When ξ 3 = 0.42 and X A = 0.25 the large particles display four well defined peaks, while the small particles have a flatter structure. Correspondingly, the velocities of the two species (lower right panel) are anti-correlated with the density profiles, the largest velocities being attained where the density displays local minima. Also the velocity profiles at such large packing fraction show an inversion, that is the larger species has the larger velocity. The other case considered for ξ 3 = 0.26, X A = 0.75 displays a less pronounced structure and nearly no oscillations in the velocity profiles.
To conclude, the density profiles are rather sensitive to the composition and show the characteristic enhancement near the walls as packing increases. On the contrary, we have not found a sensitive dependence of the density profiles on the applied load. The velocity profiles for moderate packing fractions have shapes reminiscent of the parabolic Poiseuille-like profiles, with the velocity of the larger species being smaller at low packing , in agreement with the prediction of the toy model of the previous section.

B. Coarse grained observables.
Besides considering the microscopic aspect represented by the various profiles, it is of interest to analyze some average properties, such as the volumetric flow of each component and the selectivity, quantities which are more easily measurable in experiments. We define the volumetric flow rate of each species as: where A is the area of the plates. In the upper panels of Figs. 7 and 8 we display how these quantities vary with the packing fraction, ξ 3 , for a fixed value of the mixture composition, X A . In the lower panel we display how the selectivity, S, defined as varies with the packing at fixed composition. We first observe that the volumetric flows at low densities increase almost linearly with packing . In general, the volumetric flow of the small species for equal composition (X A = 0.5) is larger for the small species, on account of the larger effective volume available. If packing increases further, the volumetric flow decreases. Such a feature has the following explanation: in the low density region, Q α increases almost linearly with packing, but as ξ 3 becomes larger it determines a quadratic increase in viscosity and thus a decrease of the velocity of the fluid in the pore. As a result, there exists an optimal value of the packing for which the volumetric flow is maximum. The peak shifts towards larger values of the packing as the pore becomes larger. Interestingly, such a non-monotonic behavior is captured also by the toy model of section IV . The results suggest that particles could be separated by their size. Hydrodynamic chromatography, instead, exploits the Poiseuille velocity profile in channels [7].

VI. CONCLUSIONS AND PERSPECTIVES
To summarize, we have proposed a new theoretical method to study mixtures under confinement by using concepts derived from both DFT and LBM. The method has been applied to recover bulk properties such as the shear viscosity, and to study inhomogeneous situations, such as the flow in a slit-like channel of molecular thickness. Whereas MD simulations provides a complete microscopic picture for molecular level mechanisms at the price of a demanding computer time, the present method offers an efficient alternative in terms of CPU-time and high numerical versatility.
We have extended the kinetic model described in refs. [34,36] to a multicomponent system. The formalism enables us to compute not only the transport coefficients of fluid mixtures, but also to investigate inhomogeneous dynamical properties. We have verified that the numerical solution is in agreement with the values of viscosity predicted by our theory. In addition, we have considered the dynamics of a mixture of particles of different size confined between two parallel plates. Under high confinement, different species display different velocity profiles. We also found non-trivial behavior when the plate separation becomes so narrow that the harshly repulsive forces among the particles and with the confining walls come into play.
We wish to comment that our approach relies on the validity of the local equilibrium approximation and is corroborated by the observation that both local thermodynamic and hydrodynamic theories appear to be reliable up to surprisingly small wavelengths. At smaller scales, both conditions break down and a transition to a genuinely molecular regime takes place. Compared to previous approaches [40][41][42], our method contains two main ideas which in our opinion render it a workable scheme: the simplified kinetic kernel (15) and the numerical scheme based on the LBM techniques. Very fast numerical solutions can be obtained along these lines and complex geometries can be handled with great simplicity and good analytical insight.
In developing the numerical algorithm via a generalized version of the Lattice Boltzmann method, we did not take into account temperature variations across the system because we wanted to maintain the methodology as simple as possible. The non isothermal extension of the algorithm is doable, but this would require to include Hermite polynomials up to fourth order (that is, 127 mesh neighbors in the collisional part of the algorithm), in order to guarantee numerical stability and physical accuracy. In the future, we plan to handle this aspect.
Before concluding we wish to add that several applications of our theory can be envisaged to the study of fluids in nanospaces. Among others we mention just a few: a) the spreading of a highly concentrated region in flow conditions and of the competition between the advection and molecular diffusion [53]; b) the transport of different molecular species in a channel of varying width and the influence of the size and/or mass on the dynamical properties [7]; c) the entropic, Asakura-Osawa interactions with walls and the interplay between advection and diffusion. When particles of different size are present one expects that the larger particles experience a larger attraction towards the walls, the so called effective entropic force [54,55]. At low Reynolds number, the scenario should not be greatly modified, but the issue is worth to be considered with some care; d) the dynamics of menisci formed in narrow slits between two coexisting liquid phases [60]; e) driven flows past chemically, physically heterogeneous substrates, a mechanism which could enhance mixing [57]; f) since the equilibrium phase diagram of substances under confinement is shifted with respect to bulk conditions, a similar dynamical effect is expected and could be used to obtain specific properties [58]; g) cylindrical pores are known to lead to strong, non-trivial size selectivity, which is highly sensitive to the pore width [59]. It would be very interesting to exploit our model in such a geometry.
Appendix A: The connection with the chemical potential In this appendix we give a proof of formula (35). To this purpose let us consider the r.h.s. of the momentum equation (10) for species α. At equilibrium the l.h.s. of eq. (10) vanishes together with the drag and the viscous contributions, so that we obtain the condition: Under the adiabatic hypothesis, we assume the non-equilibrium positional pair correlations to be given by the equilibrium pair correlations [26,27] evaluated when the instantaneous density profiles assume the values n A (r, t) and n B (r, t). The force F α,mf (r) can be expressed as where we have used the property of the hard sphere potential: in order to recast the surface integral over the sphere of radius σ αβ eq.(35) into a volume integral. Substituting the result (A1) into eq. (A2) we obtain the Born-Green-Yvon relation [25]. We now compare the r.h.s. of eq. (A1) with the exact equilibrium relation where F int is the non-ideal contribution to the free energy functional, −k B T c αβ (r, r ′ ) = δ 2 Fint δn α (r)δn β (r ′ ) is the Ornstein-Zernike direct correlation function and −∇µ α int is the force arising from the interactions with the other fluid particles on a particle located at r [43]. By equating the r.h.s. terms of eqs. (A2) and (A4, we conclude that eq. (35) holds.
Finally, notice that eq. (A4) can be rewritten as that is, in thermal and mechanical equilibrium the net force on a fluid particle at position r must vanish. This occurs when the entropic force, represented by the first term, being the resultant of the attractive and repulsive forces exerted on a molecule at r from all the other molecules (the second term on the l.h.s.) and the external force (on the r.h.s.) compensate exactly.
where k q are the nodes of the quadrature over a spherical surface [65] and W q the associated quadrature weights. We choose a 18-points quadrature that is exact for quadratic integrands. The quadrature nodes on the unit sphere are the six on-axis points (±1, 0, 0), (0, ±1, 0) and (0, 0, ±1) and twelve points along the diagonals (± √ 2, ± √ 2, 0), (± √ 2, 0, ± √ 2) and (0, ± √ 2, ± √ 2). Correspondingly, the quadrature weights are W q = 1/15 for the on-axis points and W q = 1/30 for the diagonal points. It should be borne in mind that the quadrature points do not necessarily fall on the mesh. For instance, if the the HS radius is an integer and one considers intraspecies collisions, there are only six quadrature nodes that fall on the mesh while twelve nodes along the diagonals do not. In order to sample the fluid density and velocity on off-mesh points, we employ a trilinear interpolation scheme [66]. It is clear that the accuracy of the collisional integrals improve for larger HS diameters, since the eight mesh points employed in the trilinear interpolation fall closer to the spherical surface. The same type of consideration applies for the interspecies integral and, in all circumstances, for the midpoint rule employed to estimate the radial distribution function via the Carnahan-Starling expression eq. (44). In our numerical applications, we applied a trilinear interpolation for any off-mesh quadrature nodes and midpoints for the radial distribution function, and chose HS diameters to be even numbers in order to minimize the number of off-mesh interpolations.
In solving eq. (15) with the LBM approach, we need to control the numerical error introduced by the discretization procedure over velocity, space and time variables. It is well-known that for one-component, uncorrelated fluids LBM bears three types of error as compared to the analytical solution of the incompressible Navier-Stokes equation. The first error depends on the mesh spacing ∆x and decays quadratically with ∆x. In presence of hard boundaries, and in particular with curved fluid-wall interfaces that are approximated as staircase geometries, the accuracy degrades to linear dependence with the mesh spacing. The second type of error depends quadratically on the time-step ∆t. Finally, the so-called compressibility error arises from the fact that LBM does not exactly enforce a non-zero divergence of the velocity field, with a resulting error that depends quadratically on the Mach number. As a side effect, the time-step error also degrades linearly with ∆t. These three types of error can be made arbitrarily small by changing ∆x, ∆t and the mass unit while keeping fixed the physical value of the mesh spacing, kinematic viscosity and mass density. In order to minimize the composite error, the mesh spacing is typically the parameter that is reduced, by negotiating with the numerical effort needed to resolve local morphological details and flow pattern. A concomitant effect of the space-time discretization is the effective reduction of the kinematic viscosity arising from a negative viscosity of numerical origin, − v 2 T ∆t 2 , that depends on the specific form of the lattice employed, such as D3Q19 in our case. The numerical contribution to viscosity is evaluated via a Chapman-Enskog analysis with the outcome that the effective viscosity can be controlled to arbitrary accuracy [19]. For the case of correlated dynamics, the viscosity due to the finite-time propagation enters with the same numerical value of the uncorrelated case, as shown in our previous publication [36].
The extension of the LBM to non-ideal fluids, such as the binary mixture described in this paper, carries the same types of error as for the uncorrelated case. The BGK contribution to the mixture dynamics, encoded by the kernel eq. (15), generates the same type of numerical error as for the single-component case. The issue of resolving the atomic-scale structural correlations critically depends on the rapidly varying oscillations of the density profiles, as much as on the current and higher moments profiles. The profiles need to be resolved at spatial scale smaller than the molecular radius ∆x < min α σ αα /2. When confronting with experiments performed on nanoscopic systems, the computational load needed to solve the kinetic equation requires a number of mesh points that scales as ∼ σ 3 αα . The same argument applies when resolving the temporal evolution, with a time-step that must be smaller than the typical collision time, ∆t ≤ min α σ αα /v T . Regarding the additional compressibility arising in LBM, the error plays a minor role as compared to the highly compressible nature of fluids at molecular scale. The low-Mach numbers often employed in nanofluidics, condensed-matter conditions further alleviates this problem.
In addition to the standard error terms, the computation of the collisional integrals eq. (B10) introduces a further source of error. It is a simple exercise to show that for a uniform system the diagonal points in the quadratures produce the collisional contribution to viscosity. Therefore, the numerical error on the collisional viscosity are mostly due to the off-lattice location of the quadrature points and the employed trilinear interpolation scheme. The level of accuracy of viscosity can be controlled in a systematic way and the theoretical curve can be recovered to an excellent level. Our numerical results indicate that, once the molecular spatial scale is resolved to account for the microscopic flow patterns, the error due to the quadratures becomes negligible.
Lattice nodes adjacent to the walls obey collision rules different from those characterizing bulk nodes. We have adopted the so called bounce-back collision rule which states that the velocity of a particle incident on a wall is reversed after the collision.