Challenges in modelling diffusiophoretic transport

The methodology to simulate transport phenomena in bulk systems is well-established. In contrast, there is no clear consensus about the choice of techniques to model cross-transport phenomena and phoretic transport, mainly because some of the hydrodynamic descriptions are incomplete from a thermodynamic point of view. In the present paper, we use a unified framework to describe diffusio-osmosis(phoresis), and we report non-equilibrium Molecular Dynamics (NEMD) on such systems. We explore different simulation methods to highlight some of the technical problems that arise in the calculations. For diffusiophoresis, we use two NEMD methods: boundary-driven and field-driven. Although the two methods should be equivalent in the limit of very weak gradients, we find that finite Peclet-number effects are much stronger in boundary-driven flows than in the case where we apply fictitious color forces.


Introduction
Chemical potential gradients in a bulk fluid cannot cause flow, as they do not result in net forces on any sub-volume of the fluid. The reason is that there are two ways in which the momentum of a fluid element can change: 1) due to a net externally applied force (e.g. gravity) on the particles in the volume and 2) due to a net imbalance between momentum flowing in through opposing boundaries of the volume element. But momentum flux through a boundary is what we normally call pressure. Therefore, an imbalance of momentum flux through opposing boundaries results if the pressure were not uniform.
If we consider a bulk fluid at constant pressure, and in the absence of external forces, other thermodynamic driving forces, such as gradients in T or µ, cannot cause net forces on a fluid element.
To illustrate this, consider a bulk binary system composed of N f solvent particles f and N s solute particles s. We assume that the composition is not homogeneous. Then each species is subject to a chemical potential gradient ∇µ i , for i = s, f . We consider the case that the pressure in the bulk of the fluid is constant, and for simplicity, we also assume that the temperature is constant. Although the system as a whole is not in equilibrium, we assume local thermodynamic equilibrium. We can then write the Gibbs-Duhem relation as V dP − S dT = N s dµ s + N s dµ f = 0; (1) which, at constant P and T , implies: It is often convenient to interpret a gradient in the chemical potential of species i as (minus) a force that acts on this species. The introduction of such fictitious, species-dependent "color" forces is allowed because the gradient of a chemical potential has the same effect as the gradient of a real potential acting on a given species. This is, of course, well-known for electrolyte solutions where gradients of the electrostatic potential and the chemical potential have the same effect.
Importantly, the Gibbs-Duhem equation (1) establishes a relation between the color forces: if each particle of species i is subject to a color force F i ∼ −∇µ i , then Eq. (2) expresses the fact that the net force on a fluid element vanishes.
However, contrary to what happens in the bulk, a gradient in the chemical potential of the various components in a fluid mixture can cause a net hydrodynamic flow in the presence of an interface that interacts differently with the different species in the solution. In Fig. 1, we show a flat solid wall and a binary solution composed of solutes s and solvents f . Each species interacts with the wall differently, with solutes being adsorbed preferentially at the solid surface. The adsorption creates an excess of solutes in the diffuse layer. Moreover, if there is a chemical potential gradient on the solutes ∇µ s , then they move following the thermodynamic force −∇µ s . As a result of the excess at the interface, the solute movement drives the solution flow. All this takes place within the diffuse layer, beyond which the fluid moves force-free; thus we observe the typical plug-like flow [46,28]. Such a flow, induced by chemical-potential gradients, is known as diffusio-osmosis. Other flows that are enabled by the presence of an interface are electro-osmosis and thermo-osmosis, each one having an "excess" quantity associated. The former originating from an excess of charges and the latter from the excess enthalpy at the interface. rµ s < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 7 n C 4 Z m F Q C J u n L i G D T f Z e 7 D 7 s i w = " > A A A C B 3 i c b V D L S g M x F M 3 U V 6 2 v q k s 3 w S K 4 s c x U Q Z d F N y 4 r 2 A d 0 h p L J Z N r Q J B O S j D A M / Q B / w K 3 + g T t x 6 2 f 4 A 3 6 H a T s L 2 3 r g w u G c e z m X E 0 p G t X H d b 6 e 0 t r 6 x u V X e r u z s 7 u 0 f V A + P O j p J F S Z t n L B E 9 U K k C a O C t A 0 1 j P S k I o i H j H T D 8 d 3 U 7 z 4 R p W k i H k 0 m S c D R U N C Y Y m S s 5 F / 4 Q 4 U i 6 P N 0 o A f V m l t 3 Z 4 C r x C t I D R R o D a o / f p T g l B N h M E N a 9 z 1 X m i B H y l D M y K T i p 5 p I h M d o S P q W C s S J D v L Z z x N 4 Z p U I x o m y I w y c q X 8 v c s S 1 z n h o N z k y I 7 3 s T c X / v H 5 q 4 p s g p 0 K m h g g 8 D 4 p T B k 0 C p w X A i C q C D c s s Q V h R + y v E I 6 Q Q N r a m h R Q 5 y j T F e m K L 8 Z Z r W C W d R t 2 7 r D c e r m r N 2 6 K i M j g B p + A c e O A a N M E 9 a I E 2 w E C C F / A K 3 p x n 5 9 3 5 c D 7 n q y W n u D k G C 3 C + f g E X z J o 1 < / l a t e x i t > rµ s < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 7 n C 4 Z m F Q C J u n L i G D T f Z e 7 D 7 s i w = " > A A A C B 3 i c b V D L S g M x F M 3 U V 6 2 v q k s 3 w S K 4 s c x U Q Z d F N y 4 r 2 A d 0 h p L J Z N r Q J B O S j D A M / Q B / w K 3 + g T t x 6 2 f 4 A 3 6 H a T s L 2 3 r g w u G c e z m X E 0 p G t X H d b 6 e 0 t r 6 x u V X e r u z s 7 u 0 f V A + P O j p J F S Z t n L B E 9 U K k C a O C t A 0 1 j P S k I o i H j H T D 8 d 3 U 7 z 4 R p W k i H k 0 m S c D R U N C Y Y m S s 5 F / 4 Q 4 U i 6 P N 0 o A f V m l t 3 Z 4 C r x C t I D R R o D a o / f p T g l B N h M E N a 9 z 1 X m i B H y l D M y K T i p 5 p I h M d o S P q W C s S J D v L Z z x N 4 Z p U I x o m y I w y c q X 8 v c s S 1 z n h o N z k y I 7 3 s T c X / v H 5 q 4 p s g p 0 K m h g g 8 D 4 p T B k 0 C p w X A i C q C D c s s Q V h R + y v E I 6 Q Q N r a m h R Q 5 y j T F e m K L 8 Z Z r W C W d R t 2 7 r D c e r m r N 2 6 K i M j g B p + A c e O A a N M E 9 a I E 2 w E C C F / A K 3 p x n 5 9 3 5 c D 7 n q y W n u D k G C 3 C + f g E X z J o 1 < / l a t e x i t > The preferential interaction of the solutes with a solid surface creates an excess of this species at the interface. The thermodynamic force −∇µs drives the solute motion creating a net flux due to the excess at the interface, defining the flow of the whole system. Surface-induced, "phoretic" flow phenomena are usually negligible in macroscopic channels, but can become dominant in micro or nano-scale channels, as phoretic fluxes scale as the channel diameter squared, whereas Poisseuille fluxes scale as the fourth power. From now on, we will often use the term "phoretic" transport to the wider class of surface-induced flow phenomena, even though, strictly speaking, phoresis is the phenomenon where particles move under the influence of the same gradients that can cause flow along fixed surfaces.
Simulations provide a tool to gain a better microscopic understanding of the factors that affect phoretic flows. In particular, simulations could make it possible to predict the strength of such flows based on the knowledge of the relevant intermolecular interactions. This in contrast to the more traditional descriptions that make use of hydrodynamic continuum theory and thermodynamics. Clearly, the need for quantitative understanding of phoretic transport is growing as more research focuses on nano-scale phenomena. But simulations of phoretic transport require special care, as they require approaches that differ from their bulk counterparts. Over the past years, much progress in this direction has been made. In this paper, we focus on one particular form of phoretic transport, namely diffusio-osmotic flow.
Diffusio-osmotic flow is a subject that was introduced by Derjaguin, using the language of thermodynamics and hydrodynamics. As an example, the presence of the colloid perturbs the neighbouring fluid creating a heterogeneous region close to its surface known as the diffuse layer. We consider the case that the colloid radius a is much larger than the thickness L of the diffuse layer. Derjaguin introduced this "boundary layer approximation" [12], to separate the problem into two regions: one inside and the other outside the diffuse layer. Due to the scale separation, the dynamics can be studied inside the diffusive layer. In this approximation, the diffusio-phoretic problem reduces to studying the flow of a fluid induced by a gradient of chemical potential parallel to a flat surface (see Fig. 2. a < l a t e x i t s h a 1 _ b a s e 6 4 = " u I K U d j v D l t n 2 I U k 7 6 p a a 1 5 X 6 n d 5 R E U 4 g 3 O 4 B B d u o A 4 P 0 I A W U E B 4 g V d 4 s 5 6 t d + v D + l y 0 F q x 8 5 h S W Y H 3 9 A h n Z l d E = < / l a t e x i t > L < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 t u l g Y J h G U N h 8 r j L M e m q e 1 e P 7 + E = " x y 5 g Q s w P n 6 B T h P m g 4 = < / l a t e x i t > Fig. 2. Diffusio-osmosis can be seen as diffusiophoresis under the boundary layer approximation. Rather than focusing on the movement of the colloidal particle, we focus on the fluid flow on its surface. In the case that a L, this reduces the analysis to a fluid flow on top of a flat plate, known as Derjaguin's approximation.
As was already noted by Barrat and Bocquet [5], a continuum approach is perfectly adequate to describe the hydrodynamics of fluids at a distance of more than a few molecular diameters from a surface. However, in order to estimate the magnitude of the velocity profile close to the surface, a microscopic picture is needed. In a sense, this fact is already clear from the long discussion about the meaning of the ζ-potential in electro-kinetic flows: this quantity depends sensitively on the local fluidity and molecular arrangements near a surface and typically cannot be predicted with any accuracy on the basis of macroscopic arguments alone. Moreover, even though the action of the surface is usually very local, its effect extends into the bulk as it changes the effective hydrodynamic boundary conditions [6].
Early Molecular Dynamics (MD) simulations to investigate diffusio-osmosis were performed by Ajdari & Bocquet [1], who made use of the Onsager relations (see e.g [9]) to measure the diffusio-osmotic transport coefficient indirectly by measuring the excess solute flux due to an applied pressure gradient. More recent simulations have used both equilibrium and non-equilibrium MD techniques to study diffusio-osmosis [46,27,28,31]. In addition, there have been several reports on MD simulations of diffusiophoresis of colloids [38,44] and short polymers [36].
In this article, we aim to discuss the existing theories and MD methods to study diffusio-osmosis(phoresis) in a unified framework. Our first case study is the diffusio-osmotic flow in simple planar geometry. We first derive the expression for the entropy production for the transport driven by chemical potential gradients, using non-equilibrium thermodynamics. A crucial step is to construct a consistent set of thermodynamic forces and fluxes, which allows finding the Green-Kubo expressions for the diffusio-osmotic transport coefficient. We perform simulations using Non-Equilibrium molecular dynamics applying microscopic forces that represent the effect of chemical potential gradients. Moreover, we propose an alternative route to derive a theoretical estimate for the diffusio-osmotic flow velocity. Our general expression reduces to the well-known theoretical results by Derjaguin [2] and Anderson [12] in the limit of an ideal-dilute solution in the bulk.
Later, we study colloidal diffusiophoresis. We examine a spherical particle under the influence of solutes in a binary solution. We performed simulations using two non-equilibrium techniques. We first imposed the explicit thermodynamic force driving the phoretic motion. Alternatively, we used the microscopic representation of the chemical potential gradient. We show that the hydrodynamic regime given by the Peclet number is crucially different for the two approaches.

Simulation techniques
Phoretic transport occurs in systems out of equilibrium, but if the applied driving forces are small enough, it is possible to estimate transport coefficients using linear response theory [46,31]. The main advantage of working in the linear regime is that it allows us to compute transport coefficients by studying auto or cross time-correlation in equilibrium. As was shown by Onsager [33], the transport matrix that provides the linear relation between the fluxes J i and the (thermodynamic) forces X j is symmetrical:  We show a BD-NEMD simulation box with source and sink regions where we impose, respectively, a high and a low concentration of the red particles. The FD-NEMD box can be viewed as representing the same thermodynamic state as a small region in local thermodynamic equilibrium (LTE) from the BD-NEMD simulation. In FD-NEMD, the concentration is homogeneous and a force representing the effect the chemical potential gradient is applied to each particle.

Non-Equilibrium Molecular Dynamics (NEMD)
In real systems, phoretic transport is the result of some externally imposed gradient in the thermodynamic fields (temperature, chemical or electrical potential) that determine the equilibrium properties of a system. In a simulation, one can choose to impose such inhomogeneities, but alternatively, one can apply a fictitious external field that has the same effect as these inhomogeneities. The idea behind this approach becomes clear if we consider Einstein's derivation of the relation between the diffusion coefficient D and the mobility m of a particle [13]: m = D/k B T . Einstein considered the balance between the flux of particles under the influence of a concentration gradient and the counterbalancing flux due to an external potential gradient.
The advantage of using an external field, rather than the original concentration gradient, is that the field can be kept constant in a system with periodic boundary conditions, whereas gradients due to real variations in the concentration must be periodic, to be compatible with the boundary conditions. Furthermore, we assume that different species can be driven by different fields. Fields that act specifically on particles of a given type only, are usually called color forces. Although the results of the calculations should not depend on which approach is chosen, we shall see later that, in the non-linear regimes, interesting differences appear.
In what follows, we shall refer to simulations where the fluxes are due to periodically repeated concentration differences, as Boundary-Driven Non-equilibrium MD (BD-NEMD). If the fluxes are due to constant color fields, we will refer to Field-Driven Non-Equilibrium MD (FD-NEMD).
We note that the gradients could also be created between two large but finite reservoirs, in which case periodic boundary conditions would not be needed, provided the conditions in the reservoirs are being kept fixed. However, without periodic boundary conditions, net flow along the direction of the gradients is not possible.
In BD-NEMD we create two spatially separated reservoirs in the simulation box: typically, two slabs separated by half the box-length in the x-direction. The particle concentrations in these two slabs are fixed at different values, to maintain a concentration gradient. In FD-NEMD we apply an external color force on each particle, which mimics the influence of the thermodynamic force. In Fig. 3, we show the connection between BD-NEMD and FD-NEMD simulations of bulk diffusion.

Boundary-Driven Non-equilibrium Molecular Dynamics
The most intuitive way of imposing a chemical potential gradient in a simulation is to explicitly create two reservoirs in the simulations separated by a transport region as shown in Fig. 3. In this case, the concentrations at the boundary of the transport region, define the flux within, thus the name "boundary-driven". The first simulations of systems experiencing chemical potential gradients in the context of diffusion were developed almost simultaneously by Heffelfinger and van Swol [21] and MacElroy [29]. The former authors called the method Dual Control Volume Grand Canonical Molecular Dynamics (DCV-GCMD), as it consists of two grand canonical MC (GCMC) control volumes or reservoirs embedded in an MD-NVT simulation box. The GCMC serves to keep the desired concentration in the reservoirs. The molecules flow between the two control volumes, from the source at high concentration to the sink at a lower concentration. Replenishing the particles in the reservoirs at the right rate generates a steady-state flux of particles. This step is critical, as it may give incorrect results if the MC/MD frequency is not large enough [21,4,8]. The tuning depends on the size of the reservoirs, the distance between them and the number of GCMC insertion/deletion attempted per MC step.
BD-NEDM is inherently inhomogeneous. The approach is perfectly suited to simulate microscopically inhomogeneous systems, such as the flow through nanoscopic films [40]. However, in other cases, the method has many disadvantages. As discussed before, it is difficult to tune the parameters to set up the initial concentration profile. Moreover, the use of GCMC implies that the velocity of the inserted particles must be known a priori and the method becomes problematic for fluid mixtures with large size ratio [40]. As we will discuss below, the magnitude of the gradient can lead to simulations occurring outside the linear response regime [4]. Finally, the simulations tend to be time-consuming, as they must explicitly include the reservoirs, and there is an overhead associated with the MC insertion/deletions, or at the very least swaps of particle identities.

Field-Driven Non-Equilibrium Molecular Dynamics
Simulations using FD-NEMD require the introduction of an external field mimicking the effect of a thermodynamic force. In general, this synthetic force has no clear physical interpretation, but its mechanical nature facilitates the simulation [14] .
The FD-NEMD approach has been applied extensively by [15] et al., with the body force coupling to particle variables such as the mass or the charge [43]. In the case of diffusion, Maginn et al. [30] performed NEMD using a color field, in which particles are assigned color charges according to their chemical identity. In this way, they replaced the chemical potential gradient by a color force of equal magnitude but opposite sign.
In practice, as the external field is non-conservative (i.e. as it is not the gradient of a potential), it will result in a constant dissipation in the system at steady state. Hence, the use of color fields must be combined with the use of a thermostat. In what follows, we will make use of the Nosé-Hoover thermostat, as it is a global thermostat, and it conserves linear momentum [22].
In 2001, Arya et al. [4] wrote about the use of color forces: "this method has not been widely used, perhaps because the equivalence of such a homogeneous external forcing function that drives diffusion and an actual chemical potential gradient has not been formally demonstrated". However, subsequently, Yoshida et al. [46] justified replacing the imposed gradients with a constant color-force field on the basis of linear response theory, from which it also follows that the Onsager reciprocity relations hold for phoretic transport. Han et al. [19] used a different method to simulate thermo-osmosis, which assumes that the forces on fluid elements can be computed from the gradient of the local, microscopic pressure tensor profile near a solid wall. However, as discussed in refs. [27,28,17], the stress route is problematic in an inhomogeneous system (e.g. close to a wall) as the definition of the microscopic stress tensor is not unique. Different definitions of the stress tensor lead to different estimates for the force and ultimately to different diffusio-osmotic flow velocities.
To summarise, the advantages that FD-NEMD offers over BD-NEMD are that it allows the simulation the effect of a constant chemical potential gradient under periodic boundary conditions. Moreover, we can use a homogeneous simulation box compatible with local thermal equilibrium. Lastly, we will show in Secs 3.2 and 4.3 that the use of FD-NEMD makes it possible to explore also (mild) non-linear effects.

Diffusio-osmosis
Before starting the discussion, it is worth pointing out that in the literature, concentration and chemical potential gradients are taken as equivalent driving forces for diffusion. As concentration gradients are not proper thermodynamic driving forces, we will not use them, even though they are related to chemical potential gradients. To illustrate the difference: in an ideal solution, the driving force is proportional to the gradient in the logarithm of the concentration rather than a gradient in the concentration. The language of chemical potential gradients is absolutely essential to take into account that not all gradients are independent, because of the Gibbs-Duhem relation. In the language of concentration gradients, this effect is less obvious and often assumed to be negligible [18] (for a discussion, see [37]).

Diffusio-osmosis and entropy generation
In our description of diffusio-osmosis, we consider a n-component fluid in contact with a solid surface, as shown in Fig. 4. Initially, the only thermodynamic forces acting on the system are the chemical potential gradients of each species i, ∇µ i . The fluid can be divided into two regions: the bulk, where the fluid can be considered homogeneous, and the vicinity of the (solid-liquid) interface, where the concentration of the different species at a distance z from the interface, c i (z), differs from its bulk value. This deviation from the bulk concentrations decays as the distance from the surface is increased. The reason why we first consider the expression for the entropy production is because it contains both the thermodynamic driving forces and the conjugate fluxes [35].     We start from the expression for the entropy production with no temperature gradient or chemical reactions [9],

Interface
Here Φ is the dissipation function, which has units of energy density per unit of time. It is proportional to the rate of entropy production σ s and represents the dissipation of energy by an irreversible process in a control volume [24]. The gradient in the chemical potential can be expressed as where c j indicates that the derivative is evaluated at constant concentration of the additional n − 1 species with j = i. Additionally, we know that with ν i being the partial molar volume of species i. Therefore, we can express the dissipation function in Eq (3) as The total volume flux in the system Q is defined as which is the average volume flow velocity in the system. We can then express the dissipation function as The expression in Eq. (8) is convenient as it separates the diffusive fluxes, which are Galilei-invariant, from the fluid flow, which is not. As we assume that in any infinitesimal volume element local equilibrium holds, we can use the Gibbs-Duhem relation, where N i is the number of particles of species n.
Defining the densities c i ≡ N i /V , we can rewrite Eq. (9) as: Eq. (10) establishes a general relation between the thermodynamic forces in the system at constant temperature. If we choose (∇µ i ) P,T in Eq. (8) as the independent driving forces then ∇P is fixed. Conversely, if we use ∇P as a driving force, then one of the (∇µ i ) P,T is linearly dependent on the others. The connection between thermodynamic forces (fluxes) avoids problems arising from treating them independently as discussed by Gupta et al. [18]. Note that the pressure that can be held constant in an experiment is the bulk pressure [9,26]. If we impose a bulk pressure gradient, there will be fluid flow. However, even when the pressure in the bulk of the fluid is constant, the presence of chemical potential gradients can still cause a pressure gradients at an interface. If we hold the pressure in the bulk constant (∇P = 0), we can express the solvent chemical potential (∇µ f ) P,T using Gibbs-Duhem in the bulk as, The dissipation function depends on n − 1 chemical-potential gradients, plus the term due to an explicit pressure gradient: In what follows, we focus on a two-component system, with solvent f and solute s. The dissipation function then becomes where we have defined the excess flux of solute as Finally, we can write the transport matrix connecting the fluxes with the thermodynamic forces, By including the factor 1/T in the thermodynamic forces, we can cast the entropy production in a simple bi-linear form in fluxes and thermodynamic forces. Such form is needed to derive the Onsager reciprocity relations for the transport coefficients M αβ . In practice, the factor 1/T is often absorbed in the transport coefficients.

Transport coefficients
To compute the transport coefficients M αδ in Eq. (15) using FD-NEMD, we need to represent the thermodynamic forces as fictitious mechanical forces incorporated in the Hamiltonian of the system and that can act on the particles in the fluid. Here, we recapitulate the derivation by Yoshidaet al. [46] to show that such an approach provides the mechanical route to Onsager's symmetry relations. We consider a system with N interacting particles satisfying Hamiltonian equations of motion: where F i is the force exerted on particle i by all the other particles, and F ext is the mechanical equivalent of the thermodynamic force.
For the diffusio-osmotic case, we represent all chemical potential gradients by equivalent forces F µ i on every particle of species i. To satisfy the condition of mechanical equilibrium in the bulk, the force F µ s on the solute particles must be balanced by a force F µ f on the solvent particles, such that: where N B and N B s are the total number of particles and the number of solute particles in the bulk. Eq. (18) is the mechanical equivalent of the Gibbs-Duhem equation.
Expressing everything in terms of the external force on the solutes The Hamiltonian coupling of the particles to the external driving forces is It is worth pointing out that all the sums in Eq. (20) are in the bulk B. Next, we consider a system confined in a slit. The total volume of the fluid Ω includes an interfacial region. The previous expression is still valid, giving rise to the diffusio-osmotic flow, as now there is a non-vanishing contribution from the externally applied forces F ext From linear response theory [20], we can compute the response of a given observable B to an external perturbation of the form ∆H = A(x i )F 0 = H ext as Focusing on the non-diagonal terms of the transport matrix on Eq. (15), when a chemical potential gradient is applied, the observable we want to measure is the total flux of the particles Q It is convenient to write the variable that couples to the external field aṡ Finally, using Eq. (22) we can express the total volume flux as Hence, using transport equations in Eq.(15), we can establish the connection between the thermodynamic force and its microscopic counterpart as Eq. (25) is general (i.e. it is valid for arbitrary forces). However, the Green-Kubo expression in Eq. (22) is only valid in the linear regime, in which case the fluxes that appear in the expression for the entropy production (Eq. (15)) are linear functions of the thermodynamic forces. Eq. (25) seems to differ from the result reported by Yoshida et al. [46], but this is only apparent: the discrepancy is due to an unfortunate definition for ∇µ s in ref. [46], which is only correct in the limit of infinite dilution. As a consequence, F µ s of ref. [46] is underestimated by a factor φ B f ≡ N B f /N B . We now focus on the off-diagonal term M JQ of the transport matrix. This coefficient expresses the dependence of the excess solute flux on the bulk pressure gradient. A pressure gradient exerts a force on a volume of fluid rather than on individual particles. As a first approximation, one might tend to connect the thermodynamic force acting on the system to the microscopic force as (see e.g [41,16,46,43,31]): It is important to realize that in confined geometries, and a fortiori in porous media, it may be problematic to work with local pressure gradients even though it is perfectly legitimate to consider the pressure difference between the reservoirs on either side of the system. The reason is that if the potential energies inside and outside the slit are different, ∇P would show δ-function spikes at the entrance and exit of the slit, whereas local thermodynamic equilibrium requires that all µ i s are continuous. If the properties of the slit are constant in the direction of the flux, the chemical potential gradients, and hence the color forces, are constant inside the slit. Of course, due to interactions with the wall, the fluid density may vary in the direction perpendicular to the wall. In that case, a constant force per particle creates different pressure gradients at different distances from the wall. This is not in contradiction with the statement that the pressures are the same everywhere inside the reservoirs, precisely because the local pressure may vary rapidly at the entrance and exit of the channel.
In what follows, we consider a small volume ω at a distance z ω from the wall. We obtain that the Hamiltonian coupling to the external force is therefore A(ω) = ω i∈s x i . The variable that couples to the external field F P is given bẏ Finally, using Eq. (22) we can express the excess solute flux as: By comparing Eq. (28) with Eq. (15) the pressure gradient that corresponds to a constant force per particle is given by: We thus conclude that the expressions for the transport coefficients in Eq. (24) and Eq. (28) are equivalent, as the correlation functions are symmetric in time. Thus, M JQ = M QJ , fulfilling Onsager's reciprocal relations.This result suggests that to obtain the cross-coefficients, in principle, we can apply pressure gradients or chemical potential gradients. In practice, the advantage of the latter is that they do not depend on the distance with the interface.

Local and global fluxes
It is instructive to look at the expression for the entropy production in a system between two reservoirs at different chemical potentials. We will assume that the temperature of the system is constant. In that case, the pressure in both reservoirs is a function of the chemical potentials and is therefore not an independent thermodynamic variable. The dissipation function for a macroscopic volume with chemical-potential profiles µ i (r), where i labels the n different species, is given by: where the j i (r) denote the diffusive fluxes of species i, and the integral is over the surface of the system. We focus on the practically important case that the system is in contact with two external reservoirs (denoted by I and II) that, individually, are at constant µ i . These reservoirs are not in direct contact with each other. In that case, the boundary conditions are completely specified by the µ I i and µ II i . The global dissipation function of the system is then given by where J i denotes the total flow of particles of species i from I to II, i.e., the total number of particle of species i that crosses either surface per unit time.
We can use Gauss's theorem to rewrite Eq. 31 as We can rewrite this as: We note that, in steady state, the divergences of all fluxes must vanish. Hence the second line of Eq. 33 vanishes and we are left with Note that adding the rotation of a vector field to the fluxes will not change this result, provided that the normal component of this rotation at the boundaries vanishes. Another way of saying the same thing is that Eq. 30 shows that adding any flux j that vanishes at the boundaries of the system (or, at least, is purely tangential to the boundaries), will not contribute to the entropy production. The above argument also holds for other fluxes, such as the heat flux, which, in contrast to the heat flow into and out of a system, are not uniquely defined.

Local Thermodynamic Equilibrium and the Derjaguin-Anderson theory for diffusio-osmosis
We consider again the system in Fig. 4. The mixture is at a constant temperature and, we assume a chemical potential gradient of species i in the x-direction. If the bulk fluid is incompressible, hence, the density and pressure equilibrate instantaneously in this region. Moreover, the rate of the spontaneous decay of chemical potential gradients over a distance scales as 2 /D i (D i denotes the diffusion coefficient of species i). As a consequence, chemical potential differences across the boundary layers equilibrate very quickly compared to the time scale of the diffusio-osmotic flow. Therefore, we can employ local thermodynamic equilibrium, assuming that the system is in equilibrium in the z-direction, even though a chemical potential gradient can be maintained along the x-direction. Hence, we can write the relation between the thermodynamic forces in the bulk from Eq.(10) as: where P xx refers to a component of the pressure tensor parallel to the surface. At the interface, the density profile c i (z) depends on z. The fact that µ = µ exc (z) + k B T ln c i (z) is constant across the diffusive boundary layer (and for a fixed x) implies that the excess chemical potential µ exc will, in general, depend on the distance z from the wall. At a point z within the diffusive boundary layer we can write, Once more, it is important to stress that mechanical forces in liquids can only be caused by body forces such as gravity or by pressure gradients [17]. The reason why chemical potential gradients near a surface cause fluid flow is that they induce a pressure gradient near a wall. It is the pressure gradient in Eq. (36) which moves the fluid. As the chemical potential µ is constant, we can relate the concentrations in the bulk (z → ∞) and close to the surface as thus, we can rewrite Eq. (36), where ∆µ i (z) exc = µ exc i (z) − µ exc i (∞) is the excess chemical potential due to the presence of the interface. We can now combine Eq. (38) with the Stokes equation to estimate the flow velocity in the x direction: Assuming a constant viscosity η,we get We note an important feature of Eq. 40: it is not expressed in terms of the pressure gradients, which are microscopically ill-defined in an inhomogeneous system (e.g near a surface). This in contrast a chemical potential gradients which, as we know, can be replaced by a uniquely defined force per particle. The other point to note is that we have assumed that the macroscopic creeping-flow approximation holds. This assumption is, in general, not correct, and in simulations we do not make this assumption. Using non-slip boundary conditions, and exploiting the fact that outside the diffuse layer, the velocity does not vary, we obtain the bulk velocity of the fluid Note that in fluid dynamics, the slip velocity is usually defined as the velocity at the interface where the boundary condition is imposed. However, in the present case, using a local continuum description, the slip velocity v B x is the fluid velocity in the bulk just outside the diffuse layer.
The Derjaguin-Anderson description of diffusio-osmosis [12,2] can be obtained as a special case of Eq. (41). For an ideal bulk solution, we have: Thus, we can write Eq. (41) as, If we now restrict the analysis to very dilute solutions of solute molecules s in a continuum liquid phase where we neglected the solvent contribution as e −β∆µ exc . In Derjaguin-Anderson theory φ(z) is the mean-field potential acting on solutes at a distance z from the solid surface. This potential does not only include the direct effect of the surface on the solutes, but accounts for the perturbation of the local liquid structure near a wall, as illustrated in Fig. 4. We note that using the excess chemical potential instead of φ has the advantage that the expression c i (z)/c B i = [e −β∆µ exc i (z) − 1] follows directly from the fact the chemical potential depends only on x and not on z. Hence we can write To quantify whether the net effect of this layering is an accumulation or depletion of particles near the surface, we use Gibbs's definition of the surface excess for particles of species i as (see [3]): Γ i is positive if there is net adsorption of particles on the wall, and negative in the case of depletion at the interface. Following Anderson, we define the so-called solute adsorption-length K i to the zeroth moment of the excess-concentration profile: K i can be interpreted as the thickness (positive or negative) of a layer of bulk solution that would contain the same net number of adsorbed or depleted particles. K i is obtained experimentally by equilibrium adsorption studies and it can be as large as 1 µm [3], or even much larger near a wetting transition. A second measure of the adsorption/depletion layer is given by ξ i . ξ i , which has the dimensions of length squared, is related to the first moment of the excess concentration: Derjaguin defined the characteristic extension of the diffuse adsorption layer as √ ξ i [11], and Anderson defined the characteristic length L * i 1 : Using the above definitions, we can rewrite the diffusio-osmotic velocity in Eq. (44) v where α is the concentration gradient of solutes in the bulk. Note that even when there is strong net adsorption of solutes (large K s ), L * s may be small, zero, or even of the opposite sign, depending on (c s (z) − c B s ). In other words, diffusio-osmotic flow is less sensitive to the excess concentration closest to the wall. This effect becomes very pronounced for thick adsorption layers, in particular near a wetting transition.

Simulations
There are two ways of imposing microscopic forces for diffusio-osmosis using FD-NEMD 5. As discussed above, we can mimic a chemical potential gradient by applying color forces on all particles [46]. Alternatively, and more in the spirit of hydrodynamics, we can start from the force per volume element [28], and then express the force per particle as where the force on the solutes F µ s is given by Eq. (25) and the force on the solvents F µ f is determined by imposing mechanical equilibrium in the bulk (see Eq. (18)). Both approaches should give the same flow profiles if the spatial binning used to measure the concentration distributions in Eq. (50) is the same as the one used to sample the velocity profiles.
In practice there is a difference: the force F µ i on solute and solvent particles is the same throughout the system. However, the net force per volume element is non-zero only close to the wall. In our simulations, we make use of this fact by imposing F µ ave (z) =0 at distances that are sufficiently far away from the wall for the density modulations to have decayed. The advantage of this approach is the following: due to the dynamic adsorption/desorption of solute and solvent particles near the wall, there will be small -in the thermodynamic limit: vanishingly small -composition fluctuations in the bulk. However, even though these fluctuations are small, their integrated effect is non-negligible: they would result in a flow in the direction opposite to the interface-induced flow. This spurious bulk flow is a finitesize effect, in the sense that, for a sufficiently large wall area, positive and negative density fluctuations will cancel. However, in a finite system, we need to suppress this spurious bulk flow explicitly. This we do by truncating the force per unit volume outside the interfacial region. In Sec. 4 we will consider a more complex geometry where we cannot easily work with the average force, and we will describe an alternative method to suppress the spurious bulk flow. I u E 7 s k d V C i 7 Z k / z j D C S U i 4 w g x J O b C t W L k Z E o p i R v K q k 0 g S I z x B I z L Q l K O Q S D e b Z s / h i V a G M I i E H q 7 g V P 1 7 k a F Q y j T 0 9 W a R V C 5 6 h f i f N 0 h U c O l m l M e J I h z P H g U J g y q C R R F w S A X B i q W a I C y o z g r x G A m E l a 5 r 7 k s 8 T i X F M q / q Z u z F H p Z J t 9 m w z x r N u / N 6 6 6 r s q A K O w D E 4 B T a 4 A C 1 w C 9 q g A z B 4 A i / g F b w Z z 8 a 7 8 W F 8 z l Z X j P L m E M z B + P o F C c S d w A = = < / l a t e x i t > Fig. 5. Different methods to impose microscopic forces for diffusio-omotic(phoretic) simulations. On the right-hand side, a force F µ i is applied on each particle depending on their species i. On the left-hand side, the average force for each particle F µ ave computed using Eq (50) is shown.
We benchmarked our simulations against the published results of ref. [46]. All simulations were perfomed using the LAMMPS software package [34]. Particles interact via a 12-6 Lennard-Jones potential (LJ) V LJ (r) = 4 LJ ij [(σ LJ ij /r) 12 − (σ LJ ij /r) 6 ] shifted and truncated at r = r cut , such that The indices i and j denote the particle types in our simulations: solutes (s), solvents (f ) and wall (w). We chose the same Lennard-Jones interaction for the particle pairs ss, sf , f f with LJ ij = 0 and σ LJ ij = σ 0 , such that the bulk solution is an ideal mixture. For convenience (Ockham's razor) we also use these same parameters for the wall-solvent interaction wf . We also assume that all particles have equal mass. The wall-solute interaction strength LJ ws and σ LJ ws were varied to control the degree of solute adsorption or depletion. For all interactions, r c = 2.5σ 0 . In what follows, we use the mass m 0 of all the particles (s,f and w) as our unit of mass and we set our unit of energy equal to ε 0 , whilst our unit of length is equal to σ 0 . All other units are expressed in terms of these basic units.
A snapshot of the system in contact with the confining wall is shown in Fig. 6. The initial dimensions of the simulation box are (17σ 0 , 17σ 0 , 35σ 0 ) with 7424 solution particles. The average concentration of solutes in the whole volume ofc s = 0.15. The box is periodic in the x and y directions. In the z direction, there is a solid wall at the bottom and a moving surface at the top, where particles undergo specular reflection. As the tangential momentum of particles remains unchanged upon reflection, this wall imposes "slip" boundary conditions. In contrast, lower surface consists of a layer of immobile solid atoms with the structure of the (100) surface of a face-centred cubic (FCC) lattice with lattice constant √ 2σ 0 . The interaction parameters of the solutes with the wall are (ε sw , σ LJ sw ) = (1.5, 1.5). We used a Nosé-Hoover thermostat [22] to fix k B T /ε 0 = 1.0 for all the simulations. To initialise the system, we performed 10 5 NVT MD steps, using a time step ∆t = 0.002τ . 5 × 10 5 steps were required to impose P σ 3 0 /ε 0 = 1.0 as described in [47,45]; this was achieved by allowing the box height to fluctuate, with the imposed pressure applied to the moving wall. During this process, we sampled the height in the z-direction. For all the subsequent simulations, the height was fixed at the average value of this fluctuating height. After equilibrating the system, we sampled the density distribution for all the species during 3 × 10 6 steps (see Fig. 7). The initial peak of the solvents near the wall is due to the fact that the wall-solute repulsion is stronger than the wall-solvent interaction. The migration of some solute particles towards the interface during the equilibration decreases their concentration in the bulk. Therefore, c B s <c. However, the effect is negligible for the system size and the relatively weak solvent adsorption, provided that we use the simulation technique described below. As the reflecting top surface is hard, there is also some layering of the fluid there (see Fig. 7). However, as there is no specific adsorption or depletion at the reflecting wall, it does not contribute to phoretic transport, and we can ignore it in our subsequent analysis. Liu et al. [28] have shown that the flow profiles obtained using Eq. (50) are in good agreement with results obtained applying an explicit chemical potential gradient. As explained above, the application of F µ ave (z) will not reproduce the correct diffusive fluxes, but that is not important in the context of this paper. One disadvantage of Eq. (50) for computations is that it requires knowledge of the equilibrium concentration profiles. In the present case, we have computed these profiles in a separate simulation, using a bin-width of 0.25 σ 0 . However, it would probably be better to use the "bin-less" method of refs. [7,10].
A conceptual disadvantage of working with F µ ave (z) rather than the forces per species, is that ∇ × F µ ave (z) need not be zero. In contrast, the rotation of the color forces vanishes.
The force distribution on the solution is shown in Fig. 8. From this figure it is clear that the net force per volume element is negligible for z > 4σ 0 . To obtain good statistics for the diffusion osmotic flows, we needed long simulations (10 8 time steps). We applied the computed force distributions and measured the velocity profiles in the fluid. Results in Fig. 9 show the diffusioosmotic velocity profile for ∇µ s = −0.125.We observe the plug-flow profile characteristic of diffusio-osmosis. At the interface, there is initially a steep increase in velocity due to the excess of solutes. Comparing with the benchmark, our results show a higher diffusio-osmotic velocity which comes from the fact that the color force used in ref [46] underestimates the effect of the thermodynamic force ∇µ s by a factor equal to the molar fraction of solvents in the bulk φ B f = N B f /N B . Notice that all the flow profiles are non-monotonic in z and exhibit a peak before settling down to the bulk velocity. This peak has also been observed in previous studies [46,28]. This overshoot can only be partially be described using Eq. (40). The remaining disagreement is not surprising, as Eq. (40)

Comparison with Theory
In order to evaluate the theoretical expressions (Eq. (40) or Eq. (44)) for the slip velocity, we need to compute the concentration distribution c i (z) of all species i as a function of the distance z from the wall, and the viscosity η of the solution. The former is obtained from EMD simulations, and relatively short runs are required as the equilibration in the z-direction is fast. In the simplest theoretical description, the viscosity is assumed to be independent of z, and equal to its bulk value: η(z) = η B . η B was obtained using the Green-Kubo expression [20]. Assuming that η is independent of z is a strong assumption, as we know that the fluid shows layering near the wall.
All relevant parameters in Sec. 3.4, such as Γ , K, L * , depend on moments of the concentration distributions. The integrals in the definition of these parameters are evaluated from the surface (z 0 = 0) to the bulk (z → ∞). However, on a microscopic scale, the location of z 0 is problematic, the more so as the particle-wall interactions are different for solvent and solute as they have different σ LJ . In Fig. 10 we show the density profile for each species close to the wall. We can define a distance d min i as the shortest distance to the wall where particles of species i can penetrate. This distance is different for solvent and solute: we obtain d min f = 0.55 and d min s = 1.15 for the solvents and solutes respectively. But the uncertainty in the location of the wall is not the only problem with the Derjaguin-Anderson theory: if the adsorption of one or more species on the wall is very strong, we should expect the local viscosity to become large and the strongly adsorbed layer will not contribute to diffusio-osmosis.
This non-uniqueness of the location of the boundary in Eqn. 47 makes a comparison between theory and simulation difficult. In fact, there are two problems: a) the location of the boundary is different for solutes and solvents, but more importantly, b) a direct simulation of pressure-driven Poisseuille flow in the channel shows that the viscosity close to the wall is clearly higher than the bulk value, resulting in a smaller slope of the, otherwise parabolic flow profile close to the wall. The latter effect can be seen in Fig. 11. For phoretic transport, which depends on a the adsorption or depletion in a microscopic surface layer, the problems with the definition of z 0 are serious.
We note that the most important parameter, L * , depends on the first moment of the concentration profile. Fig. 12 how strongly the integrand in Eq. (47) depends on the assumed value of z 0 . The results show that in the case that we consider it is almost meaningless to attempt a quantitative comparison between the microscopic simulations and the macroscopic theory. But that is not all. As we argued above, the Derjaguin-Anderson theory of diffusio-osmosis ignores the effect of the chemical potential gradient of the solvent. However, Eq. 40 allows us to go beyond the standard Derjaguin-Anderson theory by taking into account all chemical potential gradients. Importantly, we can determine the contribution from the different species to the velocity in Eq. (40). In the present case, it is straightforward to estimate the sign of the diffusio-osmotic velocity a priori, using thermodynamic arguments [3]. However, Eq. (40) can also deal with situations where there is a multi-component solution with competing interactions between the species and the surface. Fig. 13 shows the velocity profiles for different values of the gradient of the chemical potential ∇µ s . As explained above, the problem with the comparison with the Derjaguin-Anderson theory is twofold: we need to choose the location of the non-slip boundary, and we need to assume constant viscosity. In Fig. 13 we have computed the theoretical profiles assuming that the non-slip boundary is at d min f , where the flow velocity vanishes, and we have assumed that the viscosity is everywhere equal to the bulk viscosity. As the figure clearly shows, with these inputs, there are large discrepancies between theory and simulation. Of course, a better agreement between theory and the Derjaguin-Anderson theory can be achieved by changing our choice for the local viscosity and the location of the slip plane, but then we would be fitting rather than predicting. However, even using a local viscosity is not solving the problem, as the viscosity is not a local quantity, as its Fourier transform is wave-vector dependent [42,32].
We also note that the parabolic part of the pressure-driven velocity profile in Fig. 11 extrapolates to zero at a distance from the wall where the real flow velocity is non-zero.
Finally, we note that the Derjaguin-Anderson theory predicts an small overshoot in the velocity profile, which can be attributed to the layering in the fluid close to the wall.

The effect of finite Péclet numbers
In our discussion of diffusio-osmosis, we considered two simulation techniques: on the one hand an approach where explicit (periodic) concentration gradients are imposed (the "Boundary-driven" NEMD approach) and on the other hand, an approach where we imposed fictitious color fields that reproduce the effect of chemical potential gradients (the "Field-driven" NEMD approach). Although in principle, both methods are equivalent, we chose to use the FD-NEMD approach, as it is computationally more convenient. However, in some situations, there are large differences in simulations using FD and BD NEMD. To be more precise: the two methods are still equivalent in the limit where the gradients vanish, but non-linearities show up much more strongly in the BD-NEMD approach than in the FD method. In this section we discuss an example, namely colloidal diffusiophoresis, where these differences can be shown quite clearly.

Boundary-Driven Non-Equilibrium Molecular Dynamics
For our BD-NEMD, we used a double-control-volume semi-grand canonical algorithm. We use a semi-grand canonical ensemble that allows us to swap particles between the two reservoirs. As we consider again solvent and solute particles that are otherwise identical, all swap moves are accepted.
The size of the simulation box was (51.30 x 20.52 x 30.78) (in units of σ 0 ). A colloid was fixed in the centre of the simulation box (see Fig. 14) by placing a large Lennard-Jones particle with σ cs = σ cf = 3.23 σ 0 , where the subscript c denotes the colloid. The concentration gradient was created by using two reservoirs of particles. The source region at c B s = 0.6σ −3 0 and the sink at c B s = 0.15σ −3 0 . The difference in concentration between the reservoirs is equivalent to ∇µ s ∼ 0.06. The imposed concentration gradient is linear when the interaction of the colloid with solvent and solute is the same: ε cs = 1.0 (see Fig 15). Note that, rather than probing the steady-state notion of a colloid in a stationary fluid, we compute the (equivalent) steady state flow of the fluid past a fixed colloid.

Source Sink
Fig. 14. Dual control volume simulation box used for the boundary-driven non-equilibrium simulations. In both control volumes, the concentration for each particle species was fixed, with the sink and source indicating the low and high solute concentration regions respectively. The distance between the reservoirs is ∆x ss = 12x l and the length of the control volumes in the x direction is ∆ cv x = 3x l , where x l = 5 1/3 σ0.
As before, the MD simulations were carried out using LAMMPS, and with the same model for solvents (f ) and solutes (s). Moreover, σ cs = σ cf = 3.23 and ε cf = 1. The only difference is that we vary the interaction strength between colloid and solute ε cs to reproduce solute depletion (ε cs = 0.5) and attraction (ε cs = 2.5). We initialized the system with a solute/solvent ratio c B s /c B f = 1 and an average solution density in the box of c = 0.75σ −3 0 . We swapped particle identity in the reservoirs every 20 time steps, with a time step of ∆t = 0.05τ . We let the system equilibrate for 10 7 steps. By doing this, we achieve both equilibration of solutes around the colloid and the desired concentration gradient between the control volumes. The equations of motion were integrated using a velocity-Verlet algorithm, and we kept the temperature of the system at k B T / 0 = 1.0 using a Nosé-Hoover thermostat [22]. After equilibration, we ran 10 7 production steps to sample the flow profile around the colloid and the concentration distribution for each species.

BD-NEMD Results
In Fig. 15 we show the solute concentration profiles for different colloid-solute interactions ε cs . As soon as phoresis starts, i.e. for LJ cs = 0, the concentration gradient becomes non-linear due to advection. This is a consequence of the finite Péclet number, which in this case is of the order of P e ≈ vL/D, where v is the average flow velocity, L the distance between the reservoirs, and D the diffusion coefficient of both solute and solvent. As a result, the local concentration gradient at the location of the colloid decreases (see also [44]).
In fact, a simple argument shows that the concentration profile should become approximately exponential. To this end, we consider the fluxes of solvent and solute in a steady velocity field v. We ignore the fact that the colloid presents an obstacle. We can then write the flux of species i as a sum of a diffusive and a convective contribution, In steady state, the concentration profile must be of the form where D s /v = 1/k defines the characteristic length scale of the concentration profile. The coefficients α and β are determined by the boundary conditions at the source and sink regions of the system (see Fig. 14).
with ∆c B s = c source s − c sink s . We can define the Péclet number for the BD-NEMD simulations as P e BD = k∆x ss = v∆x ss /D s . In Fig. 16 we show P e BD for the different interactions ε cs . D s = 0.13 was computed from the mean-square displacement of the solutes in the bulk region. Even for the smallest non-zero phoretic flow P e BD is non-negligible, and the BD-NEMD simulations cannot be used to estimate diffusiophoresis.
For freely moving colloids, the effect of the finite Péclet number on the speeds of diffusiophoresis is well known [25,3]. However, the large effect of a finite Péclet number in simulations with the boundary-driven flow seems less known. In fact, Sharifi et al. [39] reported BD-NEMD simulations of diffusiophoresis, but in order to suppress the Péclet effect, they were forced to make the concentration profile piece-wise linear, which introduces unphysical sources and sinks in the diffusive fluxes throughout the simulation box . Péclet number P e BD for the diffusiophoretic flow with several colloid-solute interaction strengths εcs. Note that, even for the smallest non-zero phoretic flow velocities P e BD is larger than one.

Field-Driven Non-Equilibrium Molecular Dynamics
To carry out FD-NEMD simulations of the same model system, we used a simulation box (20.52 x 20.52 x 30.78) (in units of σ 0 ). The system was initialized using the same procedure as in the BD-NEMD simulations. To equilibrate, in this case, we imposed semi-grand canonical swap moves between s and f throughout the simulation box. We attempted to swap 10 4 particle identities every 10 steps for the first 10 5 steps, thereby generating an equilibrium distribution of solutes around the colloid, and an equimolar solution in the bulk. The equilibration step is crucial as our aim to carry out simulations under conditions where the composition of the bulk fluid is kept fixed, even as we varied the colloid-solute interaction ε cs [36].

FD-NEMD Results
As discussed previously, we represent the chemical potential gradients by equivalent external forces that are compatible with the periodic boundary conditions. As before, the forces are chosen so that there is no net force on the system as a whole. Hence, there is only one independent force to be chosen. In the present case, we chose to fix the force on the solutes F µ s . To facilitate comparison with the BD-NEMD simulations, we fixed this force such that it corresponds to a linear concentration gradient in the BD-NEMD case. This choice resulted in, F µ s = 0.06 0 /σ 0 . Following the discussion in Sec. 3.5 and bearing in mind the complex geometry of the present case, we applied color forces on the solvent and solute, rather than average forces on the fluid as in diffusio-osmosis (see Eq.(50)).
Having specified the force on the solutes, the force on the solvent particles F µ f follows from mechanical equilibrium in the bulk: N B s , N B f denote the number of solutes and solvent in the bulk region. Once the forces in the bulk are specified, we obtain the phoretic force on the colloid F µ c by imposing mechanical equilibrium in the whole system, N s , N f refer to the number of solutes and solvents in the whole system. This equation expresses the fact that there can be no net external force on the fluid: if there were, the system would accelerate without bound, as there are no walls or other momentum sinks in the system. Eq. (56) establishes a connection between all chemical potential gradients (or the corresponding microscopic forces), which must be balanced throughout the system as the phoretic flow cannot cause bulk flow.
In practice, we exploit Galilean invariance, and keep the position of the colloid fixed. As discussed before in the context of diffusio-osmosis (Sec. ??), there are inevitably fluctuations in the bulk concentrations due to exchanges between adsorbed and non-adsorbed particles. These variations would lead to unphysical velocity fluctuations in the bulk (unphysical, because in the thermodynamic limit this effect goes away), creating noise in the observed phoretic flow velocity. To suppress this effect, we could either adjust the composition in the bulk domain at every time step or recompute the forces on the solvents F µ f ) such that the external force on the bulk domain is always rigorously equal to zero (this also adjusts the force on the colloid F µ c ). We opted for the latter approach as particle swaps would affect the stability of the MD simulations.
In Fig. 17 we show the results obtained using the BD-NEMD and FD-NEMD. The first point to note is that the BD-NEMD simulations yield a phoretic flow velocity that is systematically lower than the value obtained from the FD-NEMD simulations. The underlying reason is that whereas the characteristic Péclet number in the BD-NEMD case is determined by the system size (P e ∼ Lv/D), the Péclet number for FD-NEMD is determined by the colloid size (P e ∼ σ c v/D, which in our case is about an order of magnitude smaller. We note that the dependence of the phoretic velocity on the strength of the interaction between colloid and solute is non-monotonic. The reason is that initially, ε cs increases the excess of solutes around the colloid, which, in turn, increases the phoretic velocity as expected in the linear regime. However, for large ε cs , the closest solutes to the colloid are tightly bound and lose their mobility. Hence, they stop contributing to the flow around the colloid (see a discussion in [36]).